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1. 


EXECUTIVE  SUMMARY 


This  report  provides  a  description  of  the  Covariance  Simulation  Program  (CSP) 
developed  as  a  tool  to  examine  cesium  clock  stability  and  to  investigate  the  utilization  of  Global 
Positioning  System  (GPS)  data  in  the  Omega  system  cesium  clock  synchronization  process.  The 
current  process  of  maintaining  Omega  system  network  synchronization,  that  is  synchronization  or 
each  Omega  Station  (OMSTA)  to  "Omega  System  Time",  requires  several  steps.  The  long  term 
synchronization  of  each  OMSTA  to  a  common  time  standard  has  historically  been  accomplished 
using  reciprocal  path  Omega  readings.  Such  readings  are  processed  at  a  central  location  to  provide 
phase  corrections  required  for  each  transmitter  signal.  Because  of  telecommunication  limitations, 
these  long  term  phase  corrections  are  ordered  by  the  central  location  at  weekly  intervals.  Short 
term  synchronization  of  the  network  is  achieved  using  the  available  cesium  beam  frequency 
standards,  or  clocks,  located  at  each  OMSTA.  A  correspondence  between  Omega  System  Time 
and  Universal  Coordinated  Time  (UTC)  can  also  be  achieved  since  4  of  the  5  northern  hemisphere 
Omega  Stations  can  obtain  Loran-C  measurements,  which  provide  a  link  to  UTC.  Using  this  data, 
the  central  location  can  regulate  the  entire  system  so  that  Omega  system  time  provides  a  close 
approximation  to  UTC. 

The  process  of  re-establishing  "coarse"  synchronization,  should  synchronization  be  lost 
due  to  a  long  term  power  outage  for  example,  is  difficult  given  the  nature  of  Omega.  To  facilitate 
the  restoration  of  "coarse"  synchronization,  GPS  receivers  were  installed  at  all  Omega  stations  in 
the  mid-1980’s.  GPS  is  currently  an  experimental  system  and  cannot  be  relied  on  as  an  exclusive 
source  of  synchronization  for  the  Omega  network.  Motivated  in  part  by  the  greater  accuracy  of 
GPS  relative  to  Loran-C,  and  GPS  world-wide  availability,  GPS  timing  data  from  one  station  was 
incorporated  into  the  synchronization  process  in  1985.  The  results  from  this  initial  experience  led 
to  an  expanded  program  involving  GPS  data  from  three  other  stations  in  1988.  The  cumulative 
results  of  these  efforts  indicate  that  GPS  could  play  a  more  prominent  role  in  the  Omega 
synchronization  process  once  GPS  becomes  fuliy  operational.  The  precise  manner  in  which  GPS 
data  is  used  must  be  carefully  assessed. 

The  CSP  was  developed  to  aid  in  this  assessment  process.  The  main  element  of  the 
CSP  is  a  six  state  Kalman  filter  designed  to  process  GPS/cesium  and  cesium/cesium  phase  offset 
measurements,  and  to  thus  provide  cesium  phase  offset  estimates.  The  measurements  that  drive 
the  filter  consist  of  field  recorded  data  obtained  from  the  Omega  Stations  in  Norway  and  Hawaii. 
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These  two  stations  provide  special  data  recording  capabilities  that  allow  the  raw  measurements  to 
be  recorded.  To  fulfill  the  role  of  an  analysis  tool,  the  CSP  allows  the  user  great  flexibility  in 
defining  system  parameters,  most  notably  the  measurement  update  interval.  Using  the  CSP,  several 
hypothesized  measurement/  corrections  incorporation  schedules  were  examined.  Results  from  the 
initial  CSP  runs,  presented  in  the  form  of  plots  both  within  this  report  and  in  the  accompanying 
addendum,  confirm  the  integrity  of  the  basic  filter  design.  The  phase  offset  residual  plots  are  key 
to  evaluating  the  impact  of  varying  measurement  incorporation  intervals  on  the  reliability  of  the 
CSP  filter  computed  phase  offset  estimates.  These  phase  offset  residuals  are  basically  the 
difference  between  the  CSP  computed  phase  offset  estimates  and  the  measured  phase  offsets.  Thus 
they  provide  a  means  of  measuring  the  "error"  in  the  CSP  estimates,  relative  to  the  GPS  data, 
between  the  scheduled  measurement  incorporations.  These  "errors"  were  observed  to  grow  only 
slightly  as  the  interval  between  measurement  updates  was  lengthened.  In  all  instances  these 
maximum  errors  remained  below  the  300  nanosecond  level. 
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INTRODUCTION 


2.1 


BACKGROUND 


The  Omega  Navigation  System  provides  a  world-wide  navigation  capability  through  the 
use  of  eight  widely  dispersed  transmitters  operating  in  the  Very  Low  Frequency  (VLF)  band.  The 
system  user  employs  phase-of-arrival  measurements  to  obtain  position  information.  The  proper 
operation  of  the  Omega  system  therefore  requires  careful  synchronization  of  all  eight  stations  to 
a  common  time  standard.  Each  transmitting  station  contains  an  ensemble  of  three  cesium  beam 
frequency  standards  which  are  used  to  control  the  short  term  synchronization  of  the  network. 
Longer  term  synchronization  has  historically  been  achieved  by  the  processing  of  daily  reciprocal 
path  Omega  readings  at  a  central  location.  Because  of  the  telecommunications  delays  involved, 
phase  corrections  associated  with  this  longer  term  control  are  ordered  by  the  central  location 
weekly.  The  effect  of  this  control  scheme  is  to  maintain  what  is  termed  "internal 
synchronization."  This  concept  describes  the  relationship  of  the  timing  of  any  individual  station  to 
Omega  system  time.  Since  4  of  the  5  northern  hemisphere  Omega  stations  can  obtain  Loran-C 
system  measurements,  and  since  Loran-C  measurements  provide  a  link  to  Univeral  Coordinated 
Time  (UTC),  the  central  location  can  also  steer  the  entire  system  so  that,  ideally,  Omega  system 
time  equals  UTC  (ignoring  the  effects  of  leap  seconds).  This  process  maintains  the  "fine” 
synchronization  of  the  system. 


The  nature  of  Omega  is  such  that  it  is  difficult  to  establish  "coarse"  synchronization  if. 
for  example,  a  long  power  interruption  causes  it  to  be  lost.  GPS  receivers  were  placed  at  all  the 
Omega  stations  in  the  mid-1980’s  to  greatly  facilitate  'he  restoration  of  coarse  synchronization, 
should  it  ever  be  lost. 


Beginning  in  1985,  GPS  timing  data  was  used  as  an  input  to  the  synchronization  process 
at  one  OMSTA,  as  an  auxiliary  input  emulating  Loran-C.  After  proving  effective  in  this  role  for 
over  2  years,  GPS  data  was  incorporated  into  the  process  at  three  other  OMSTAs.  The  results 
indicate  that  once  GPS  becomes  an  operational  system,  it  could  play  a  greater  role  in  Omega 
synchronization.  However,  the  exact  nature  of  the  role  must  be  carefully  studied.  The  initial  step 
in  this  study  is  described  herein. 
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SCOPE  AND  PURPOSE 


In  essence  the  study  described  in  this  report  examines  the  performance  to  be  expected 
if  the  Omega  network  synchronization  were  derived  exclusively  from  the  cesium  and  GPS  readings 
available  at  a  given  station.  A  program  called  the  Covariance  Simulation  Program  (CSP)  was 
developed  to  characterize  the  performance.  Chapter  3  of  the  report  provides  a  brief  description 
of  the  Omega  station  equipment,  concentrating  on  the  source  of  the  data  used  in  this  study  and  its 
relationship  to  other  equipment.  Chapter  4  describes  the  development  of  the  CSP,  the  mam 
element  of  which  is  a  6  state  Kalman  filter  designed  to  estimate  cesium  phase  offsets  based  on 
cesium/GPS  and  cesium/cesium  measurements.  It  also  provides  details  of  the  filter  algorithm  and 
a  description  of  how  CSP  can  be  used  to  produce  simulated  results  under  new  synchronization 
methods.  An  overview  of  the  CSP-user  interface  is  also  provided.  Additional  details  of  the  CSP 
system  operation  are  contained  in  a  User's  Guide  (Reference  1).  Chapter  5  presents  and  discusses 
the  results  from  initial  testing  of  the  CSP  uang  OMSTA  Norway  and  Hawaii  data.  The  purposes 
of  these  tests  were  to  confirm  the  integrity  of  the  basic  filter  design  and  to  study  the  behavior  of 
the  cesium  clocks  located  at  these  stations.  A  complete  set  of  data  plots  is  provided  in  an 
addendum  to  this  report  (Reference  2). 


DISCUSSION  OF  GPS/OMEGA  TIMING  DATA 


3.1  REVIEW  OF  RELEVANT  OMEGA  STATION  EQUIPMENT 

This  section  provides  a  brief  description  of  those  pieces  of  the  Omega  system  relevant 
to  this  study.  Presently,  all  transmitting  stations  are  equipped  with  DA'HJM  Model  9390-5007  GPS 
Time/Frequency  monitors.  In  their  current  capacity,  the  GPS  measurements  are  used  as  pan  of 
a  process  for  synchronization  of  the  Omega  system  epoch  to  UTC.  The  eight  transmitting  stations 
observe,  record,  and  repon  phase  difference  measurements  of  UTC  minus  online  Omega  Signal 
Format  Generator  (OMSFOG).  The  individual  station’s  internal  delay  constant  is  then  applied 
to  obtain  UTC  minus  station  epoch. 

3.2  INDIVIDUAL  TRANSMITTING  STATION  INPUTS 

Omega  time  at  each  station  is  derived  from  three  Hewlett-Packard  CAQI-5061A 
Cesium  Beam  Frequency  Standards,  which  are  part  of  Unit  2  of  the  AN/FRN-30  Timing-Control 
Set.  Each  frequency  standard  drives  an  OMSFOG.  The  block  diagram  in  Figure  3.2-1  depicts  the 
configuration  of  equipment  for  measurement  of  Omega  epoch  with  the  GPS  implementation.  The 
three  frequency  standards  and  their  associated  OMSFOGs  are  designated  A,  B,  and  C.  as  shown 
in  Table  3.2-1  for  OMSTAs  Norway  and  Hawaii.  One  frequency  standard  (and  its  associated 


Figure  3.2-1  Configuration  of  Equipment  for  Measurement 
of  Omega  Epoch  by  GPS 
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TABLE  3.2.1 

FREQUENCY  STANDARD  -  OMSFOG  CONFIGURATION 


STATION 

FREQUENCY  STANDARD 

OMSFOG 

FUNCTION 

CESIUM 

DESIGNATOR 

UNIT 

UNIT 

DESIGNATOR 

Online 

486 

A 

2A1 

3A2 

A 

Norway 

Secondary 

017 

B 

2A2 

3A3 

B 

Primary 

124 

C 

2A5 

3A4 

C 

Online 

529 

A 

2A1 

3A2 

A 

Hawaii 

Primary 

554 

B 

2A2 

3A3 

B 

Secondary 

349 

C 

2A5 

3A4 

C 

OMSFOG)  is  designated  the  online  unit.  The  two  remaining  frequency  standards  are  designated 
primary  and  secondary  and  are  held  in  reserve  in  case  of  failure  of  the  online  unit.  Table  3.2-1 
shows  the  primary  and  secondary  cesium  designations  at  OMSTA  Norway  are  reversed  from  those 
at  OMSTA  Hawaii.  This  system  setup  will  affect  the  order  of  measurements  incorporated  into  the 
CSP  as  discussed  in  Section  4.3  of  the  report. 

It  is  important  to  reiterate  that  the  current  synchronization  process  uses  reciprocal  path 
Omega.  Loran-C  and  GPS  readings  for  long  term  synchronization.  In  this  study  we  will  concentrate 
only  on  GPS  and  cesium  readings  taken  at  each  station.  The  study  is  facilitated  by  the  fact  that 
special  equipment,  known  as  the  United  States  Naval  Observatory  Remote  Data  Acquisition  System 
(USNO  RDAS).  has  been  installed  at  the  OMSTAs  Norway  and  Hawaii.  The  USNO  RDAS  block 
diagrams  are  shown  in  Figures  3.2-2  and  3.2-3.  At  these  stations,  the  relative  phase  of  all  three 
frequency  standards  is  continuously  monitored  and  recorded  by  an  HP9915B  computer.  An 
HP59307  Very  High  Frequency  (VHF)  switch  is  used  to  compare  1  pulse  per  second  (PPS)  signals 
from  all  three  cesiums  on  channels  A 2,  B1  and  B2  to  the  1  PPS  signal  from  the  DATUM  GPS 
receiver  on  channel  Al.  Every  morning  a  USNO  operator  connects  to  a  Hewlett-Packard  computer 
through  a  3451  modem  and  collects  the  Precise  Time/Time  Interval  (PTTI)  data  saved  thus  far  in 
the  current  week.  Once  a  week  the  data  banks  of  Hewlett-Packard  computers  at  OMSTAs  Norway 
and  Hawaii  are  cleared  and  the  data  recording  cycle  starts  anew.  At  USNO,  a  Hewlett-Packard 
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1000  computer  is  used  to  save  the  PTTI  data  for  collection  periods  exceeding  one  week.  Selected 
portions  of  this  data,  discussed  in  the  next  section,  were  utilized  in  the  CSP. 


3.3  DATA  FORMAT 

The  USNO  RDAS  was  used  to  obtain  the  input  data  for  the  CSP.  Figure  3.3-1  is  an 
actual  sample  of  the  PTTI  data  collected  by  USNO.  Each  data  record  field  depicted  in  Figure  3.3-1 
for  OMSTA  Norway  is  described  in  Table  3.3-1.  Table  3.3-2  contains  similar  record  field 
descriptions  for  PTTI  data  provided  by  OMSTA  Hawaii.  Of  20  available  data  fields  the  following 
six  (appearing  with  asterisks  in  Tables  3.3-1,  3.3-2,  and  Figure  3.3-1)  data  types  are  utilized  by  the 
CSP  directly: 

1.  Modified  Julian  Day  (MJD) 

2.  GMT  Time  of  the  Day  (TOD) 

3.  Time  difference  between  GPS  time  and  secondary  cesium  clock 
(GPS-CS017)  for  OMSTA  Norway. 

Time  difference  between  GPS  time  and  primary  cesium  clock  (GPS-CS554)  for 
OMSTA  Hawaii. 

4.  Time  difference  between  GPS  time  and  primary  cesium  clock 
(GPS-CS124)  for  OMSTA  Norway. 

Time  difference  between  GPS  time  and  secondary  cesium  clock  (GPS-CS349)  for 
OMSTA  Hawaii. 

5.  Time  difference  between  online  and  secondary  cesium  clocks 
(CS486-CS017)  for  OMSTA  Norway. 

Time  difference  between  online  and  primary  cesium  clocks  (CS529-CS554)  for 
OMSTA  Hawaii. 

6.  Time  difference  between  online  and  primary  cesium  clocks 
(CS486-CS124)  for  OMSTA  Norway. 

Time  difference  between  online  and  secondary  cesium  clocks  (CS529-CS349)  for 
OMSTA  Hawaii. 
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Data  Acquisition  System  (OMSTA  Norway) 


Figure  3.2-3  USNO  Remote  Data  Acquisition  System  (OMSTA  Hawaii) 
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Figure  33-1  A  Sample  of  Pi’l  l  data  collected  by  USNO  RDA  System 

for  Omega  Station  Norway 


*  These  record  fields  are  the  direct  inputs  to  the  CSP 


3-6 


Table  33-1 

Covariance  Simulation  RDAS  Timing  Measurements 
(OMSTA  Norway) 


Record  Field 

No. 

Station  Configuration 
Channel  Description 

Description  of  the  Data  Field 

1 

NA 

Modified  Julian  Date  (MJD)  of  the 
data  being  reported  in  this  record* 

2 

NA 

GMT  Time  of  the  Day* 

3 

Al-Bl 

GPS-CS017* 

4 

A1-B2 

GPS-CS124* 

5 

A1-B3 

GPS-LORl 

6 

A1-B4 

GPS-LOR2 

7 

A2-B1 

CS486-CS017* 

8 

A2-B2 

CS486-CS124* 

9 

A2-B4 

CS486-LOR1 

A2-B4 

CS486-LOR2 

11 

A3-B1 

OMSFOG-CS017 

12 

A3-B2 

OMSFOG-CS124 

13 

A3-B3 

OMSFOG-LOR1 

14 

A3-B4 

OMSFOG-LOR2 

15 

NA 

CLOCK-CS017 

16 

NA 

CLOCK-CS124 

17 

NA 

CLOCK-LOR1 

18 

NA 

CLOCK-LOR2 

19 

NA 

CS486-LORM 

20 

NA 

CS486-LORX 

_ ; 

*These  record  fields  are  the  direct  inputs  to  the  CSP. 
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Table  33-2 

Covariance  Simulation  RDAS  Timing  Measurements 
(OMSTA  Hawaii) 


Record  Field 

No. 

Station  Configuration 
Channel  Description 

Description  of  the  Data  Field 

1 

NA 

Modified  Julian  Date  (MJD)  of  the 
data  being  reported  in  this  record* 

2 

NA 

GMT  Time  of  the  Day* 

3 

Al-Bl 

GPS-CS554* 

4 

A1-B2 

GPS-CS349* 

5 

A1-B3 

GPS-LOR1 

6 

A1-B4 

GPS-LOR2 

7 

A2-B1 

CS529-CS554* 

8 

A2-B2 

CS529-CS349* 

9 

A2-B4 

CS529-LOR1 

10 

A2-B4 

CS529-LOR2 

11 

A3-B1 

OMSFOG-CS554 

12 

A3-B2 

OMSFOG-CS349 

13 

A3-B3 

OMSFOG-LOR1 

14 

A3-B4 

OMSFOG-LOR2 

15 

NA 

CLOCK-CS554 

16 

NA 

CLOCK-CS349 

17 

NA 

CLOCK-LOR1 

18 

NA 

CLOCK-LOR2 

19 

NA 

CS529-LORM 

20 

NA 

CS529-LORX 

These  record  fields  are  the  direct  inputs  to  the  CSP. 
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It  should  be  noted  that  it  is  necessary  to  pre-process  some  of  the  data  before  it  is  used  as 
a  measurement  in  the  CSP.  Julian  day  and  GMT  time  of  day  are  extracted  from  the  MJD  field  in 
the  following  manner: 

47475  .2083333 

47475  =  Modified  Julian  Day 
=  10  November  1988 

.2083333  =  GMT  Time  of  the  Day 
=  .2083333  x  24 
=  0500  Hours 

Additionally,  it  is  necessary  to  convert  the  raw  timing  data  to  microseconds  and  decode  it. 
The  following  algorithm  is  used  to  perform  this  conversion: 

•  If  a  measurement  is  less  than  0.5  seconds,  it  is  simply  multiplied  by  106  to 
obtain  a  phase  (time)  measurement  in  microseconds 

•  If  a  measurement  is  greater  than  or  equal  0.5  seconds,  then  first  1  is 
subtracted  from  it,  and  then  it  is  multiplied  by  106. 

For  the  specific  timing  data  records  examined  within  this  report,  each  of  the  records  were 
reviewed  to  identify  any  anomalies  (e.g.,  data  dropouts)  due  to  the  RDAS  equipment  problems. 
These  anomalies  were  edited  out  of  records  to  create  records  suitable  for  CSP  processing. 

A  final  note  regarding  Tables  3.3-1  and  3.3-2  is  in  order.  Note  that  though  OMSFOG  data 
is  available,  it  is  not  used  in  this  study.  The  OMSFOG  data  contains  the  effects  of  corrections 
inserted  by  the  current  synchronization  process.  By  using  only  the  measurements  obtained  before 
the  OMSFOGs,  we  are  able  to  produce  simulations  which  are  independent  of  the  current  control 
process. 


4. 


COVARIANCE  SIMULATION  PROGRAM 


This  section  describes  the  Kalman  filter  developed  to  investigate  cesium  clock  stability.  The 
filter  is  designed  to  provide  cesium  phase  and  frequency  offset  estimates  based  upon  GPS  timing 
measurements.  Section  4.1  discusses  the  models  for  the  two  primary  physical  elements  in  the 
model:  the  cesium  clocks  and  the  GPS  receiver.  Section  4.2  develops  the  filter  model  and  Sections 
4.3  and  4.4  discuss  the  various  run  options  provided  by  the  program. 

The  filter  described  here  may  be  used  as  a  tool  to  investigate  hypothesized  synchronization 
schemes  that  might  incorporate  GPS  data.  Theoretical  questions  that  may  be  examined  include 
questions  such  as  how  frequently  must  GPS  data  be  processed  to  keep  the  cesium  phase  error  at 
an  acceptable  level.  The  CSP  was  developed  to  provide  the  user  the  flexibility  to  examine  different 
measurement  sets  and  measurement  update  schedules.  To  set  the  stage  for  the  filter  discussion, 
a  brief  overview  of  key  features  of  the  cesium  synchronization  process  is  provided  below. 

Cesium  frequency  standards  have  excellent  short  term  stability,  but  gradually  drift  so  that 
after  extended  intervals  their  phase  outputs  will  exceed  tolerance  limits.  To  maintain  tolerance 
limits,  some  form  of  clock  corrections  must  be  used.  Under  the  current  control  method,  new  system 
phase  corrections  for  each  station  are  derived  and  applied  once  a  week.  The  major  hindrance  to 
more  frequent  corrections  is  the  telecommunications  burden  associated  with  the  collection  of 
reciprocal  path  Omega  readings.  If  we  were  to  rely  exclusively  on  GPS  and  cesium  readings,  more 
timely  corrections  might  be  feasible  since  all  the  measurements  needed  to  derive  corrections  for  a 
given  station  would  be  obtained  at  that  station.  Though  numerous  other  implementation/quality 
assurance  factors  would  have  to  be  considered,  we  can  explore  some  theoretical  aspects  of  the 
problem. 

4.1  CLOCK  MODELS 

4.1.1  Cesium  Clocks 

Cesium  clocks  are  characterized  by  excellent  short  term  stability.  Once  synchronized,  such 
clocks  are  able  to  provide  highly  accurate  time  outputs  for  extended  periods.  Gradually,  however. 
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errors  in  the  time  outputs  increase  to  unacceptable  levels;  frequency  errors  integrate  into  phase 
errors  and  a  slow  growth  in  the  phase  offset  errors  appears.  At  such  times  some  form  of 
resync*1  ronization,  using  new  measurement  data,  is  required. 

To  develop  an  estimation  algorithm,  some  model  of  the  growth  of  such  errors  must  be 
considered.  The  model  described  here  uses  two  clock  states:  phase  offset  and  frequency  offset, 
representing  a  bias  and  a  drift  in  the  cesium  clock,  respectively.  Frequency  offsets  propagate  in 
time  into  phase  offsets,  and  random  effects  occurring  at  both  the  frequency  and  phase  levels  also 
corrupt  the  phase  offset.  The  description  provided  here  is  based  upon  Reference  3. 

The  phase  or  time  offset  of  a  transmitter  relative  to  UTC  (in  tisec)  is  given  at  time  tn+1  by* 
[&#]»+!  =  [«<«„  +  iTlsfl„  +  [»*]„  (4.1-1) 

and 

[Sf]n+1  =  [Sf]„  +  [w,]„  (4.1-2) 

where 

8<b  =  cesium  phase  offset  in  nsec, 

<5f  =  cesium  frequency  offset  in  ^sec/day, 

w^  =  zero-mean  white  noise  sequence  corrupting  the  time  offset 
(5$)  caused  by  uncorrelated  fluctuations  in  cesium  clock 
frequency, 

wf=  zero-mean  white  noise  sequence  corrupting  the  time  drift  rate 

(<5f)  caused  by  uncorrelated  fluctuations  in  cesium  clock 
frequency-rate, 

AT  =  discrete  time  step  (1/2  day). 

A  block  diagram  representation  of  these  equations  is  shown  in  Figure  4.1-1. 

*The  subscript  n  denotes  the  discrete  time  at  which  the  phase  offset  is  computed.  Thus  [<5d>]n  = 
5</>(tn)  where  tn(n=  1,2,3,...)  are  the  computational  times  separated  by  the  discrete  time  step  AT. 


Figure  4.1-1  Covariance  Simulation  Cesium  Clock  Error  Model 

This  model  shows  that  any  initial  time  drift  rate  (cesium  frequency  offset)  will  produce  a 
corresponding  time  offset  that  grows  linearly  with  time.  Any  random  fluctuation  in  the  cesium 
clock  frequency,  w^,  results  in  a  random  walk  timing  error  that  grows  in  proportion  to  the  square 
root  of  time.  The  RMS  value  of  w^  is  taken  in  this  study  to  be 

aw6  =  0.0189  /xsec  (or  crW(i)  =  3.6  x  10‘4  nsec2)  . 

The  white  sequence  wf  is  included  in  the  clock  error  model  to  roughly  approximate  the  long  term 
random-walk-type  behavior  in  the  nominal  time  drift  rate  (cesium  frequency  offset)  observed  in 
empirical  data.  Each  discrete  jump  in  the  cesium  frequency  offset  of  magnitude  A5f  occurring  after 
an  N  day  interval  is  approximated  by  a  series  of  small  jumps  of  magnitude  A5f/2N  occurring  twice 
a  day.  The  RMS  value  of  wf  appropriate  for  this  approximation  is  given  by 

<7wf  =  A5f//ZN  . 
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For  a  worst-case  jump  of  0.3  /x sec/day  and  a  worst-case  time  interval  of  15  days,  the  above  equation 

yields 

¥ 

a =  0.0 55  nsec/ day  (or  a2^  =  0.3  x  10'2  (/xsec/day)2). 

The  values  for  both  aW(^  and  a ^  presented  in  this  section  are  consistent  with  Reference  3. 

4.1.2  GPS  Receiver/GPS  Time  Outputs 

The  receiver  employed  at  the  Omega  sites  is  a  DATUM  Model  9390  GPS  Time/Frequency 
monitor.  In  the  time  and  frequency  monitoring  mode  of  operation,  the  receiver  outputs  time  pulses 
synchronized  to  UTC,  and  measures  time  and  frequency  inputs  from  an  external  time  and  frequency 
source.  In  contrast  to  the  short  term  stability  provided  by  stand-alone  cesiums,  GPS  receiver  time 
outputs  are  characterized  by  excellent  long-term  stability.  GPS  receiver  time  outputs  are  derived 
via  filtering  of  pseudorange  measurements,  with  the  resulting  time  accuracy  dependent  on  several 
error  sources.  GPS  time  errors  are  typically  modeled  by  biases  and  additive  white  noise. 

In  the  context  of  this  effort,  a  model  of  GPS  time  output  errors  is  required  to  develop  a 
model  of  the  GPS/cesium  difference  measurements.  Specifically,  the  uncertainty  in  the  GPS/cesium 
measurements  will  be  a  function  of  the  uncertainty  in  the  GPS  time  data  and  the  inherent 
uncertainties  in  the  measurement  process  itself.  In  this  study,  the  error  in  the  GPS  measurements 
is  considered  as  a  white  noise  sequence  with  standard  deviation 

c 7q  =  0.05  (/xsec)  (or  ac2  =  0.0025  jxsec2). 

This  value  has  been  selected  on  the  basis  of  empirical  data. 

Note  that  no  attempt  has  been  made  to  model  biases  in  the  GPS  time  outputs.  The  errors 
are  assumed  to  be  zero  mean  white  noise  characterized  by  the  variance  given  above.  This  is  not 
a  severe  limitation.  Owing  to  the  long  term  stability  of  GPS  time,  infrequent  corrections  of  the 
outputs  might  keep  any  biases  at  acceptable  levels.  The  consideration  of  techniques  for  removing 
biases  from  the  GPS  timing  outputs  is  beyond  the  scope  of  the  current  effort. 


4.2  KALMAN  FILTER  ALGORITHM 


This  section  provides  a  discussion  of  a  Kalman  filter  developed  to  estimate  time  errors  in 
cesium  clocks  using  GPS-based  time  measurements.  A  basic  discussion  of  the  algorithm  is  provided 
in  Section  4.2.1.  Sections  4.2.2  and  4.2.3  describe  the  user  provided  CSP  inputs  to  the  model.  The 
two  input  parameter  files  required  to  run  the  CSP,  PARM.IN  and  KFTABL.IN,  are  covered  in 
Sections  4.2.2  and  4.2.3,  respectively. 

4.2.1  Filter  Algorithm 

A  six  state  Kalman  filter  was  developed  to  provide  estimates  of  the  phase  and  frequency 
offsets  for  each  of  the  three  cesiums  using  available  measurement  data.  The  specification  of  the 
filter  is  straightforward  once  the  system  dynamics  and  observation  model  have  been  defined.  A 
discussion  of  the  system  dynamics  and  the  observation  model  used  here  are  presented  in  Sections 
4.2. 1.1  and  4.2.1. 2,  respectively.  For  convenience,  the  standard  set  of  Kalman  filter  equations 
(Reference  4)  are  summarized  in  Table  4.2-1.  The  two  major  steps  in  this  process,  the 
propagation  step  and  the  measurement  update  step,  are  discussed  in  subsection  4.2.1.3  and  4.2.1. 4. 
respectively. 


TABLE  4.2-1 

KALMAN  FILTER  ALGORITHM 


! 

i 


I 

| 

| 

I 

! 


I 


Optimal  Estimate 

x  =  4>  x  +n.j  Propagation 

*  +n  =  *  'n  +  K*  izn'Hn  *  '„)  UPdate 

K„  =  P'„HTn  [H„P  „HTn  +  RJ'1  Weighting 

Estimation  Error  Covariance 

P'„  =  <t,P+n.i‘*>T  +  Q„.,  Propagation 

P+n  =  [I-KnH„]  P  „[l-KnH„]T  +  K„R„K„T  Update 


I 


A 

The  main  objective  of  the  filter  is  to  provide  an  estimate,  denoted  5  0  ,  of  the  actual 
phase  offset  50.  The  state  estimation  error  x  n  at  time  tn  is  defined  as  the  difference  between  the 

A 

optimal  estimate  x  n  and  the  actual  system  state  xn,  i.e, 
x  n  =  x  n  ‘  *11  • 

Letting  E  denote  the  expected  value  operator,  the  error  covariance  matrix  P  is  defined  as 
Pn  =  E(  x  n  x  nT]  • 

Note  when  n=0  this  represents  the  initial  covariance. 

The  remaining  matrices  appearing  in  the  algorithm  are  process  noise  defined  as 
Q  =  E  [wn  wTn] 
and  measurement  noise  defined  as 
R  -  E  [V„  vT„). 

These  matrices  model  the  various  error  processes  that  affect  system  synchronization  and  must  be 
specified  prior  to  estimating  the  system  state  vector  Xj,.  Q  and  P0  depend  on  the  calibration  and 
performance  characteristics  of  the  cesium  frequency  standards.  The  measurement  covariance 
matrix  R  depends  on  the  magnitude  of  the  uncorrelated  measurement  errors. 

4.2. 1.1  System  Dynamics 

The  system  states  are  defined  to  be  the  phase  and  frequency  offsets  for  each  of  the 
three  cesiums:  online,  primary  and  secondary.  The  system  state  is  written  as 

x  =  [50lt  Sfj,  S02,  Sf2,  503,  5f3]T 
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where  the  indices  1,  2  and  3  correspond  to  the  online,  secondary  and  primary  cesiums  for  OMSTA 
Norway.  For  OMSTA  Hawaii,  the  index  2  corresponds  to  the  primary  and  the  index  3  to  the 
secondary.  This  convention  corresponds  to  the  order  of  the  measurements  in  the  data  file.  For 
clarity,  the  remainder  of  the  section  follows  the  Norway  convention. 


The  state  propagation  equation  has  the  form 


xn  =  +  wn.j 


The  form  of  the  state  transition  matrix  follows  from  the  discussion  in  Section  4.1.1,  and  is  given  by 


<J>= 


A  0  O' 
0  A  0 
0  0  A 


where  A  is  a  two-by-two  matrix  of  the  form 


1  AT 
0  1 


The  driving  white  noise  sequence  w  also  follows  from  Section  4.1.1, 
w  =  [w*lf  w^.  wd>2,  w*2,  w^3,  wf3]T  . 

4.2. 1.2  Observation  Model 

The  observation  model  relates  the  system  states  to  the  measurement  values.  In  the  current 
application  there  are  two  general  classes  of  measurements  with  a  total  of  five  different 
measurement  types.  The  GPS/cesium  phase  difference  measurements  are  of  three  types: 
GPS/online,  GPS/secondary  and  GPS/primary.  The  cesium/cesium  phase  difference  measurements 
are  of  two  types:  online/secondary  and  online/primary.  The  model  for  one  type  of  each  class  will 
be  discussed  here;  the  remaining  types  follow  by  analogy. 
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It  is  important  to  point  out  that  the  displayed  value  of  phase  offset  shown  in  this  report,  in 
an  absolute  sense,  has  no  particular  relationship  to  the  time  error  of  the  current  Omega  system. 
Synchronization  errors  such  as  shown  in  Figure  4.2-1  and  other  figures  in  this  report  are  obtained 
from  the  readings  of  the  cesium  outputs,  before  the  OMSFOG  corrections  are  applied.  Thus,  the 
values  of  the  displayed  phase  measurements  represent  the  "coarse"  frequency  offset  of  the  cesium, 
propagated  by  the  time  since  the  USNO  RDAS  equipment  initialization.  The  other  contributors 
to  the  displayed  phase  offsets  can  be  the  equipment  resets  or  failures  with  the  resultant  timing 
jumps.  If  a  decision  were  made  to  build  a  synchronization  system  critically  dependent  upon  the 
RDAS  equipment,  special  operational  procedures  and  redundant  equipment  would  be  employed 
at  the  OMSTAs.  The  goal  of  this  analysis  is  to  make  an  initial  estimate  of  the  system  complexity 
and  to  provide  a  method  to  determine  the  required  measurements/correction/incorporation  schedule 
to  achieve  a  required  performance  level.  Therefore,  for  the  purposes  of  this  study,  the  behavior 
of  the  USNO  RDAS  measurements,  their  variation  from  time  to  time,  is  of  the  most  importance. 
The  absolute  value  of  the  displayed  timing  data  has  no  real  significance. 

The  general  linear  relation  between  the  state  and  the  measurement  is  given  by, 

z  =  Hx  +  v. 

In  this  equation,  z  is  the  actual  measurement,  H  is  the  observation  vector  relating  the  state  and  the 
measurement,  and  v  is  the  observation  error,  assumed  to  be  uncorrelated.  The  covariance  of  v  is 
called  the  measurement  noise  matrix  (the  R  matrix). 

For  the  GPS/cesium  difference  measurements,  in  particular  the  GPS/online  cesium 
measurements,  the  model  is  as  follows. 

Z  =  -  5<t>  +  VGC  ' 

where  z  is  the  measured  GPS/online  cesium  phase  difference,  and  86  is  the  phase  offset  in  the 
online  cesium.  Comparing  with  the  general  linear  relation  given  above  and  state  definition,  it  is 
evident  that 


H  =  [-1.  0,  0,  0,  0,  0]  . 
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The  term  vGC  represents  the  combined  error  in  the  GPS  receiver  time  and  the  error  inherent  in 
the  measurement  process,  that  is, 

VGC  =  VG  +  VGCM  * 

where 

vG  =  error  in  GPS  time  outputs 

VGCM  =  error  induced  by  the  GPS/cesium  measurement  process. 

The  error  vGCM  may  be  further  specified  as  the  sum  of  the  error  in  measuring  the  GPS  receiver 
time  outputs  (vGM)  and  the  error  in  measuring  the  online  cesium  time  outputs  (vci),  such  that 

VGC  =  VG  +  VGM  +  VC1  • 

Assuming  all  errors  are  uncorrelated,  it  follows  that 

^GC  =  °^G  +  °^GM  +  ^Cl  • 

Based  on  experimental  data  and  experience,  we  can  select  representative  values  for  these  variances. 
The  measurement  noise  variance  for  the  GPS  time  estimates  is  set  at  0.0025  jisec2  (or  a  standard 
deviation  of  0.05  /xsec).  The  value  for  can  be  assumed  to  be  negligible  in  comparison  to  cTq. 
The  value  for  <rci  is  also  negligible  with  respect  to  ct^q,  but  is  included  here  since  it  is  relevant  for 
the  cesium/cesium  measurement  processes.  Under  these  assumptions,  aGC  takes  the  form 

a2 Gc  =  °^g  +  ^ci  ~  °^G  =  0-0025  /xsec2  ' 

The  reasonableness  of  this  value  (crGC  =  0.05  ^isec)  may  be  qualitatively  verified  by  examining 
Figure  4.2-1.  This  figure  shows  a  plot  of  the  GPS/cesium  measurement  for  the  online  cesium  at 
OMSTA  Norway  for  a  9  day  period.  The  assumed  standard  0.05  ^.sec  deviation  for  the 
measurements  is  shown  via  the  dashed  lines  on  the  plot.  Better  than  55%  of  the  data  points  fall 
within  these  bounds.  The  assumed  variance  is  a  reasonable  starting  point  for  this  model. 


For  the  cesium/cesium  difference  measurements,  specifically  the  online/primary  cesium,  the 
model  is  as  follows: 

z  =  +  50j  -  50 3  +  vci3 

where  z  is  the  measured  online/primary  cesium  phase  difference,  50  j  and  503  are  the  online  and 
primary  phase  offsets,  respectively,  and  vci3  is  the  error  induced  by  the  measurement  process.  The 
subscripts  1  and  3  reference  the  online  and  primary  clocks,  respectively.  The  statistics  of  vci3 
depend  on  several  factors  including  characteristics  of  the  hardware  devices  involved  in  the 
measurement  process.  Lacking  precise  information,  the  error  in  the  difference  measurement  is 
assumed  to  be  the  sum  of  the  errors  in  the  measurement  of  each  clock.  Specifically, 

VC13  =  VC1  +  VC3 

and,  assuming  the  errors  are  uncorrelated, 

^€13  =  ^Cl  +  °^C3  ■ 

4.2. 1.3  Time  Propagation 

In  between  measurements  the  state  estimates  and  the  covariance  matrix  propagate  as 
specified  by  the  system  dynamics.  The  state  estimate  propagates  as 


^(tk)  =  <D5?(tk.1), 


or,  more  specifically  for  our  case, 


5  0(tk)  =  50(tk_i)  +  AT  5  f  (tjj.j)  ,  and 


Sf(tk)  =  5f(tk.,). 
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Figure  4.2-1  Measurement  Phase  Offsets  (Norway) 


The  error  covariance  must  also  be  propagated  forward  in  time  according  to  the  dynamics 
equations.  In  general, 

P-(tk)  =  ®  P+(tk.,)  KT  +  Q(tk.,). 

The  growth  in  the  error  covariance  is  seen  to  be  comprised  of  a  deterministic  part,  $P+<I>T,  and 
a  term  accounting  for  random  disturbance,  the  process  noise  matrix  Q.  The  process  noise  matrix 
Q  follows  from  the  cesium  clock  model.  Q  is  a  diagonal  matrix  with  diagonal  elements,  given  by: 

Qn  =  Q33  =  Q55  =  > 

Q22  =  ^44  =  ^66  =  ^wf  • 

4.2. 1.4  Measurements  and  Measurement  Update 

When  measurement  data  becomes  available  it  may  be  used  to  update  the  state  estimates 
and  covariance  matrix.  The  Kalman  filter  processing  algorithm  provides  the  means  to  incorporate 
these  measurements  in  a  statistically  optimal  sense.  The  updated  state  estimate  is  equal  to  the 
a  priori  estimate  plus  the  product  of  the  Kalman  gain  and  the  measurement  residual,  i.e., 

x+  =  x  *  +  K(z-Hx') 

where. 


K  =  Kalman  gain, 

H  =  observation  vector, 
z  =  measurement, 

/v 

z  -  H  x  *  =  predicted  measurement  residual, 

^  , 

x  =  measurement  updated  estimate  (a  posteriori  estimate),  and 

A 

x  *  =  propagated  estimate  (a  priori  estimate). 


The  optimal  Kalman  gain  is  given  by 
K  =  PHt(HPHt  +  Ry1 
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where 


R  =  measurement  noise  matrix, 

P  =  covariance  matrix,  and 

HPHt  +  R  =  residual  variance. 

The  uncertainty  in  the  estimate,  expressed  in  P,  is  reduced  by  the  incorporation  of  the  new 
measurement  data 

P+  =  P‘  -  KHP\ 

Within  the  code,  this  measurement  update  equation  is  implemented  in  the  algebraically  equivalent, 
but  numerically  more  stable,  form 

P+  =  (I  -  KH)  P-(I  -  KH)t  +  KRKt. 

In  this  section  we  have  discussed  a  specific  Kalman  filter  which  defines  the  estimation 
process  required  for  this  study.  In  the  following  sections  we  will  detail  the  CSP  software  file 
structure  which  allows  the  user  to  define  the  parameters  of  the  estimation  process.  The  user  must 
define  the  measurement  and  process  noise  parameters  and  the  initial  state  and  covariance  estimate 
for  each  desired  estimate.  These  parameters  are  defined  in  a  file  called  PARM.IN,  which  is 
discussed  in  Section  4.2.2  below.  The  user  must  also  define  the  measurement  set  to  be  used,  the 
time  interval  between  measurements  and  define  the  filter  execution  and  output  process.  These  are 
defined  in  a  file  called  KFTABL.IN,  which  is  discussed  in  Section  4.2.3. 

4.2.2  Defining  System  Parameters  (PARM.IN) 

In  this  section  a  discussion  of  the  specific  filter  parameters  is  provided.  In  addition,  the 
specific  parameter  values  used  in  the  test  runs  are  provided.  The  parameters  are  related  to  the 
CSP  parameter  input  file  PARM.IN;  through  this  file  the  program  user  may  accept  default  values 
or  adjust  the  parameters  from  their  predefined  values.  The  User’s  Guide  (Reference  1)  provides 
a  detailed  discussion  of  this  file. 
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The  measurement  noise  and  process  noise  parameters,  along  with  the  initial  state  and 
covariance  estimates,  will  determine  the  filter’s  behavior.  The  process  noise  and  measurement  noise 
are  defined  on  the  basis  of  general  physical  considerations.  A  general  description  of  the  PARM.IN 
file  and  a  sample  PARM.IN  file  used  for  the  Norway  simulation  runs  are  shown  in  Table  4.2-2. 
The  various  elements  to  this  table  are  discussed  in  the  following  subsections. 

4.2.2. 1  Measurement  Noise  Matrix  (R) 

The  measurement  error  variances  for  the  five  types  of  difference  measurements  being 
considered  are  derived  from  assumed  values  for  the  each  element.  For  example,  if  the 
measurement  is  the  time  difference  between  the  GPS  and  primary  cesium  clock  times  then  the 
measurement  error  variance  value  is  determined  as 

cr2  =  +  o2^  . 

Appropriate  measurement  noise  variances  for  the  other  measurement  types  are  formed  in  a  similar 
fashion.  The  assumed  error  variances  for  the  different  timing  devices  are: 

<?-Q  =  0.0025  Msec2  GPS, 

crci  =  0.000049  /x sec2  online  cesium, 

cr2 Q2  -  0.0001  /xsec2  secondary  cesium,  and 

cr2^  =  0.000064  /x sec2  primary  cesium. 

The  value  of  o^q  has  been  established  as  discussed  in  Section  4.2.I.2.  Values  for  o^q  have  been 
chosen  to  represent  a  range  of  uncertainties  that  can  be  reasonably  expected,  but  more  importantly, 
to  demonstrate  filter  performance  under  varying  conditions.  The  parameters  of  measurement  noise 
variances  are  set  by  the  program  user  through  the  PARM.IN  file  parameters  SIG.  Specifically, 

’SIG  =  ’,  SIG  GPS,  SIG  ONL,  SIG  SEC,  SIG  PRI 
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where. 


SIG  GPS  =  standard  deviation  in  GPS  time  measurement  errors  (^seO 
SIG  ONL  =  standard  deviation  in  online  cesium  measurement  errors  (yxsec), 

SIG  PRI  =  standard  deviation  in  primary  cesium  measurement  errors  (ftsec),  and 
SIG  SEC  =  standard  deviation  in  secondary  cesium  measurement  errors  (ftsec). 


The  values  shown  in  Table  4.2-2  for  the  OMSTA  Norway  simulation  correspond  to  the  numerical 
values  detailed  above. 


TABLE  4.2-2 

PARAMETER  FILE  CONTENTS 


RECORD 

NAME 

VALUE  (UNITS) 

1 

1 

CLK 

Initial  time  offset  estimate  (fisec),  initial  frequency 
estimate  (/xsec/day). 

2 

TIME 

Time  interval  parameter  used  in  calculation  of  Q  matrix 
elements  (days). 

3 

COV 

Mean  square  uncertainty  in  phase  offset  estimate  (/isec)2, 
mean  square  uncertainty  in  frequency  offset  estimate 
(jjl  sec/day)2. 

4 

ERR 

Largest  discrete  jump  in  cesium  frequency  offset 
(/isec/day),  corresponding  measurement  time  interval  (days). 

5 

Q  COV 

Mean  square  frequency  fluctuation  (jj. sec)2,  mean  square 
frequency-rate  variance  (/tsec/day)2. 

6 

SIG 

Root  mean  square  time  measurement  error  Qu.sec).  GPS 
or  cesium  clocks 

’CLK 

=  ’,0.3570,  0.011461,  -0.7380.  0.119371,  1.21,  -0.050453 

TIME 

=  ’,0.5,  0.5,  0.5 

’COV 

=  ’,0.0001,  0.0016,  0.000225,  0.0016,  0.0000,  0.0016 

’ERR 

=  ’,0.03,  15.0,  0.03,  15.0,  0.03,  15.0 

’QCOV 

=  ’,0.00036,  0.0,  0.00036,  0.0,  0.00036,  0.0 

’SIG 

=  ’,0.050,  0.007,  0.010,  0.008 
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4.2.2.2  Process  Noise  Matrix  (Q) 


The  process  noise  matrix  Q  is  intended  to  reflect  uncertainties  in  the  form  of  the  system 
propagation  equations  (e.g.  unmodeied  states,  injection  of  random  noise).  The  form  of  the  Q 
matrix  follows  from  the  discussion  of  the  cesium  clock  physics  briefly  considered  in  Section  4.1.1. 
The  Q  matrix  is  defined  to  be  diagonal  with  the  diagonal  elements  Qll,  Q33,  Q55  corresponding 
to  o-2^,  and  Q22,  Q44,  Q66  corresponding  to  ct2^  .  The  program  user  may  adjust  these  process 
noise  tuning  parameters  through  the  P ARM.  IN  file  parameter  Q  CON. 

Q  CO V  =  ’Qll,  Q22,  Q33,  Q44,  Q55,  Q66’ 

where, 

Qll  =  online  clock  frequency  noise  variance,  and 

Q22  =  online  clock  frequency  rate  noise  variance. 

The  pairs  Q33,  Q44  and  Q55,  Q 66  are  similarly  defined  for  the  secondary  and  primary  clocks 
respectively. 

The  user  may  also  specify  the  Q22,  Q44,  Q66  terms  through  the  use  of  the  PARM.IN  file 
parameter  ERR.  This  allows  those  elements  to  be  defined  according  to  the  equation  for  a ^ 
provided  in  Section  4.4.1,  that  is 

crwf  =  ASf  /VN/TlME 

The  elements  of  ERR  are 

ERR  =  ’ERR1,  ERR2,  ERR3,  ERR4,  ERR5,  ERR6’ 

where, 

ERR1  =  ASf  for  online  cesium,  and 
ERR2  =  N  days  for  online  cesium. 

The  TIME  parameter  is  used  to  designate  the  interval  between  discrete  jumps  in  the  frequency 
offsets.  Following  Reference  3,  it  is  assumed  such  jumps  occur  "twice  a  day".  The  nominal  value 
for  TIME  is,  thus,  1/2  day.  The  remaining  parameter  ERR3,  ERR4  and  ERR5,  ERR6  are  similarly 
defined  for  the  secondary  and  primary  clocks.  (NOTE:  If  this  method  is  used  to  define  Q22.  Q44 
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and  Q66,  then  these  terms  must  be  set  to  zero  in  Q  COV,  that  is  Q  COV  must  have  the  form  Q 
COV  =  ’Qll  ,  0  ,  Q33,  0,  Q55,  0  ’). 

4.2.23  Initial  Error  Covariance  Matrix  (PJ 

In  typical  filter  applications,  the  initial  state  of  a  system  is  known  only  roughly;  a  fairly  high 
degree  of  uncertainty  may  exist.  The  diagonal  elements  of  the  initial  covariance  matrix  (PQ)  should 
be  chosen  to  reflect  this  uncertainty.  The  physics  of  the  problem,  as  well  as  experimentation,  lead 
the  filter  designer  to  choose  appropriate  values.  The  covariance  simulation  tool  allows  the  user  to 
select  these  initial  variances.  Ideally,  the  uncertainties  selected  should  reflect  the  uncertainty  in  the 
initial  conditions.  Large  variances  should  be  selected  if  there  is  a  high  degree  of  uncertainty  in  the 
initial  estimate  to  help  ensure  that  the  filter  will  converge  rapidly.  Small  variances  may  be  used  if 
the  initial  state  is  known  with  a  high  degree  of  certainty,  to  avoid  overshoots  and  oscillations. 
Selecting  variances  that  are  too  small  will  increase  the  time  required  for  the  filter  to  converge. 
Typically,  a  conservative  filter  design  will  use  large  variances  as  default  values  but  allow  the 
operator  to  override  these  defaults  if  initial  data  is  known  to  be  highly  accurate. 

The  program  user  may  select  these  values  through  the  PARM.IN  file  parameter  COV. 
Specifically, 

COV  =  ’PI  1,  P22,  P33,  P44,  P55,  P66’ 

where. 

Pll  =  Initial  uncertainty  (variance)  in  online  clock  phase  offset  estimate 
P22  =  Initial  uncertainty  (variance)  in  online  clock  frequency  offset  estimate 

the  pairs  P33,  P44  and  P55,  P66  are  similarly  defined  for  the  secondary  and  primary  clocks. 

For  the  test  cases  presented  in  Section  5,  the  diagonal  elements  of  the  initial  error 
covariance  matrix  were  set  to  100  ns2  (online  cesium),  0  ns2  (primary  cesium),  and  150  ns2 
(secondary  cesium).  The  initial  variances  corresponding  to  the  frequency  offset  states  were  set  to 
1600  (ns/day)2  for  each  of  the  three  clocks.  This  particular  selection  of  values  is  somewhat  arbitrary 
and  is  intended  in  part  to  examine  the  behavior  of  the  filter. 
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4.2.2.4  Initial  State  Estimate  (x  J 

Since  one  of  the  primary  purposes  of  this  study  is  to  examine  the  steady  state  behavior  of 
the  filter,  and  the  available  data  records  are  relatively  short,  it  is  desirable  to  begin  with  initial  state 
estimates  approximating  the  actual  state.  For  purposes  of  testing  the  performance  of  the  filter,  the 
initial  filter  estimates  were  derived  by  direct  processing  of  the  data.  That  is,  measurement  data  was 
examined  offline  to  determine  suitable  initial  estimates  for  the  phase  offset  and  frequency  offset 
for  each  of  the  three  clocks.  The  best  performance  is  achieved  if  the  initial  estimate  and  the  initial 
covariance  matrix  are  consistent. 

The  user  may  adjust  the  initial  state  estimates  through  the  PARM.IN  file  parameter  CLK. 
Specifically, 

CLK  =  ’PHI,  FR1,  PH2,  FR2,  PH3,  FR3’ 
where, 

PHI  =  Initial  phase  offset  estimate  for  online  clock, 

FR1  =  Initial  frequency  offset  estimate  for  online  clock. 

Likewise  the  pairs  PH2,  FR2  and  PH3.  FR3  are  similarly  defined  for  the  secondary  and  primary 
clocks. 


4.2.3  Defining  System  Parameters  (KFTABL.IN) 

The  second  major  input  file  involved  in  the  CSP  is  called  "KFTABL.IN".  This  file  is  used 
to  control  the  simulation  process.  It  determines  which  measurement  set  is  used,  the  time  interval 
between  measurements,  and  the  output  and  termination  functions.  There  are  four  types  of  events 
used  by  the  CSP:  update  (UPDT),  propagate  (PROP),  USER,  and  STOP.  The  operator  specifies 
parameters  for  each  of  these  events.  During  the  update  event  the  covariance  matrix  is  updated  to 
incorporate  sensor  measurements.  The  propagate  event  allows  the  covariance  to  propagate  through 
a  time  interval.  The  USER  event  exits  the  Kalman  filter  program  and  enters  the  output  routines. 

Finally,  the  STOP  event  terminates  the  simulation.  A  sample  event  file  is  shown  in  Figure  4.2-2. 
This  file  shows  three  GPS/cesium  measurements  processed  during  a  31-day  time  period.  The  filter 
process  propagates  estimates  on  a  half-day  basis,  performing  an  update  every  three  days. 
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PROP(CNST,STAT)  FROM  0.0  TO  31.0  BY  0.5 
UPDT( VARY, GPS  1,STAT)  FROM  0.0  TO  31.0  BY  3.0 
UPDT(VARY,GPS2,STAT)  FROM  0.0  TO  31.0  BY  3.0 
UPDT(VARY,GPS3,STAT)  FROM  0.0  TO  31.0  BY  3.0 
USER(LTUS)  FROM  0.0  TO  31.0  BY  0.5 
STOP  AT  31.0 _ 

Figure  4.2-2  Sample  Event  File  (KFTABL.IN) 


Events  are  scheduled  by  order  of  time  sequence  labeling  and  by  order  of  appearance  in  the 
event  file.  In  other  words,  if  several  events  are  scheduled,  they  will  be  executed  in  the  order  as  they 
appear  in  the  event  file.  The  sole  exception  is  the  PROP  event  from  tn  to  tn+j  which  is  defined 
at  tn.  It  is  executed  after  all  other  events  at  tn  and  before  any  other  event  at  tn+1. 

The  UPDT  event  identifies  the  measurement  set  used  in  the  CSP.  The  first  argument  in 
the  update  call  is  always  VARY,  indicating  a  varying  measurement  matrix.  The  second  argument 
defines  the  measurement  set.  Any  single  update  can  consist  of  five  measurements  as  shown  in 
Table  4.2-3.  The  last  argument  is  always  STAT,  forcing  the  Kalman  filter  to  propagate  the 
estimated  phase  and  frequency  offset  errors  to  the  next  sample  time  point. 

The  USER  event  is  a  general  purpose  exit  from  the  CSP  to  perform  auxiliary  computations 
on  the  covariance  analysis  in  progress.  In  the  CSP  it  is  used  to  generate  a  LOTUS  output  file. 
This  plot  file  is  a  list  of  records,  each  one  containing  time,  phase  estimate,  frequency  estimate  and 
RMS  uncertainties  in  the  estimates  A  more  thorough  discussion  of  the  filter  ouputs  and  LOTUS 
output  file  format  is  contained  in  Sections  4.3  and  4.4. 
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TABLE  4.2-3 

RDAS  MEASUREMENT  SET 


UPDT 

MEASUREMENT  DESCRIPTION 

GPS1 

Difference  between  GPS  and  on-line  cesium 

GPS2 

Difference  between  GPS  and  secondary*  cesium 

GPS3 

Difference  between  GPS  and  primary*  cesium 

CLK1 

Difference  between  on-line  and  secondary*  cesiums 

CLK2 

Difference  between  on-line  and  primary*  cesiums 

The  STOP  event  provides  the  program  the  ability  to  terminate  normally  without  hitting  the 
end  of  the  file  marker.  No  records  after  the  first  STOP  are  executed,  permitting  the  user  to  modify 
or  debug  a  simulation  easily. 


4.3  SIMULATION  OUTPUTS 

This  section  describes  the  outputs  provided  by  the  CSP.  The  reader  is  referred  to  the  User's 
Guide  (Reference  1)  for  additional  details  on  the  procedures  required  to  generate  these  outputs. 
The  primary  outputs  provided  by  the  program  are: 

•  State  Estimates  :  estimated  phase  and  frequency  offsets 

•  Filter  Indicated  RMS  Errors  (i.e.,  square  root  of  filter  variances) 

•  Phase  Offset  "Residuals" 


*  Primary  and  secondary  cesium  measurements  are  reversed  in  Hawaii  RDAS  timing  data. 
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The  main  outputs  of  the  filter  are  the  phase  and  frequency  offset  estimates.  Each  time  a 
measurement  is  provided  to  the  filter  these  values  are  updated  to  reflect  the  new  data.  In  the 
various  scenarios  examined  in  Chapter  5  and  the  addendum,  the  interval  at  which  measurements 
are  provided  to  the  filter  is  varied  between  0.5  days  and  7  days.  Between  these  updates  the 
estimates  are  maintained  by  propagation  according  to  the  state  transition  matrix.  For  the  Norway 
scenarios  the  interval  between  propagation  steps  is  0.5  days,  while  for  the  Hawaii  scenarios  this 
interval  is  1.0  day  (the  reason  for  this  difference  is  discussed  in  Section  5.2). 

The  diagonal  elements  of  the  error  covariance  matrix  P  provide  an  indication  of  the  filter’s 
determination  of  the  error  in  the  provided  estimate.  Small  variances  indicate  the  filter  is  confident 
of  the  state  estimate.  For  convenience,  the  filter  indicated  RMS  errors  (i.e.  the  square  root  of  the 
variances)  are  provided  as  outputs.  The  filter  indicated  RMS  errors  for  a  particular  phase  offset 
estimate  will  be  at  a  minimum  just  after  a  measurement  has  been  incorporated.  Between 
measurements  this  uncertainty  will  grow  for  two  reasons.  The  uncertainties  in  the  frequency 
estimate  will  propagate  into  uncertainties  in  the  phase  offset  estimates.  The  uncertainty  will  also 
increase  due  to  the  presence  of  random  disturbances  modelled  in  the  process  noise  matrix  Q. 

One  of  the  major  aims  of  good  filter  design  is  to  ensure  that  filter-indicated  RMS  errors  are 
consistent  with  the  actual  errors  in  the  state  estimates.  Ideally,  one  would  like  "truth"  data  during 
the  filter's  testing  period  to  verify  performance.  In  this  case  we  would  hope  that  the  filter's 
estimates  of  the  states  (phase  and  frequency  offsets),  and  of  the  RMS  errors  in  these  states,  would 
be  consistent.  That  is,  the  actual  error  in  the  state  estimate  should  be  roughly  equivalent  to  the 
corresponding  diagonal  element  of  the  error  covariance  matrix.  As  with  most  filter  applications  the 
"truth"  data  required  to  compute  the  actual  error  in  the  state  estimate  is  not  available.  Therefore, 
to  gain  some  insight  into  the  performance  of  the  filter,  phase  offset  "residuals"  are  also  provided 
as  an  output  by  the  CSP.  The  phase  offset  residual  for  a  clock  is  defined  here  as  the  difference 
between  the  measured  GPS/cesium  phase  difference  and  the  best  available  CSP  estimate  of  that 
cesium  phase  offset.  The  best  available  estimates  between  measurement  updates  are  the 
propagated  phase  estimates.  Further,  at  a  measurement  update  time,  two  estimates  are  available: 
an  a  priori  estimate,  propagated  forward  from  the  previous  cycle,  and  an  a  posteriori,  measurement 
updated  estimate.  This  latter  quantity  reflects  the  incorporation  of  the  new  measurement  data  and 
is  therefore  a  better  estimate  than  the  a  priori  estimate.  In  forming  the  residuals  at  the 
measurement  update  times,  the  a  posteriori,  measurement  updated,  estimate  will  be  used.  The 
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phase  offset  residuals  defined  here  are  provided  for  analysis  purposes  only  and  are  not  used  by  the 
filter  in  determining  the  estimates.  The  CSP  computes  these  phase  offset  residuals  at  each 
propagation  step  (set  at  0.5  days  for  Norway  and  1.0  days  for  Hawaii). 


This  definition  is  provided  to  avoid  confusion  with  the  term  "predicted  residual"  commonly 

A 

used  in  filter  applications.  The  term  predicted  residual  refers  to  the  quantity  z  -  H  x  *  defined  in 
Section  4.2. 1.4  and  represents  the  difference  between  an  actual  measurement  used  by  the  filter  and 
the  filter’s  predicted  measurement  based  on  the  a  priori  state  estimate.  The  predicted  residual  is 
defined  only  at  the  measurement  update  steps.  Thus,  if  the  measurement  update  interval  is  5  days, 
the  predicted  residual  is  computed  every  5  days.  Note  also  in  the  case  where  the  measurement 

/A 

used  by  the  filter  is  a  cesium/cesium  difference  measurement,  the  predicted  residual,  z  -  H  x  * ,  is 
a  qualitatively  different  quantity  than  the  phase  offset  residual  defined  in  the  previous  paragraph. 
The  predicted  residuals  will  not  be  shown  in  this  report,  instead,  the  "phase  offset  residual”  as 
defined  above  will  be  selected. 

4.4  LOTUS  OUTPUT  FILE 

This  section  discusses  the  actual  file  containing  the  various  filter  outputs.  A  single  run  of 
a  CSP  results  in  a  generation  of  LOTUS  output  file.  At  the  beginning  of  the  filter  execution,  the 
user  is  asked  to  specify  the  file  name  that  will  be  associated  with  the  output  LOTUS  file.  A  brief 
discussion  of  the  LOTUS  file  is  provided  to  familiarize  the  user  with  the  filter  applications. 

The  LOTUS  output  file  is  used  to  generate  a  set  of  simulation  plots  that  describe  the  state 
of  cesium  phase  and  frequency  offsets,  plus  the  accuracy  of  the  estimation  process.  Table  4.4-1 
defines  the  data  elements  in  the  output  file  for  OMSTA  Norway.  Figure  4.4-1  contains  an  example 
of  a  LOTUS  worksheet  generated  from  one  of  the  LOTUS  output  files. 

The  data  is  divided  into  three  major  categories.  Worksheet  columns  A  and  B  identify  the 
Julian  day  and  day  counter  from  the  beginning  of  the  RDAS  timing  data  collection  process.  The 
value  of  the  Julian  day  changes  only  with  next  measurement  update  event.  Thus,  an  example  shown 
in  Figure  4.4-1  can  be  recognized  as  a  3-day  update  rate  configuration  run.  In  other  words,  the 
measurements  get  incorporated  into  CSP  at  3  day  intervals:  0,  3,  6,  9,  etc.  Columns  C  through  H 
contain  CSP  state  estimates:  online  clock  phase  offset  estimate,  online  clock  frequency  offset 
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estimates  and  offline  clocks  respective  state  estimates.  Note  that  the  order  of  the  clock  data  (i.e., 
online,  secondary  and  primary)  is  controlled  by  the  format  sequence  of  input  (NORWAY.DAT) 
file  (refer  to  Section  3.3  of  this  report).  Columns  I  through  N  contain  uncertainties  (la)  of  the  state 
variables  contained  in  columns  C  through  H.  The  order  of  the  RMS  values  corresponds  to  the 
order  of  the  state  variables.  Finally,  columns  O  through  S  contain  measurement  residuals. 
Columns  containing  residuals  of  measurements  not  used  in  the  present  configuration  run  are  set 
to  zero. 
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TABLE  4.4-1 

LOTUS  OUTPUT  DATA  FILE  CONTENTS 
(NORWAY) 


A 

Data  Point  Modified  Julian  Date 

B 

Data  Point  Time 

C 

Phase  Offset  -  Online  Cesium  Clock 

D 

Frequency  Offset  -  Online  Cesium  Clock 

E 

Phase  Offset  -  Primary  Cesium 

F 

Frequency  Offset  -  Primary  Cesium 

G 

Phase  Offset  -  Secondary  Cesium 

H 

Frequency  Offset  -  Secondary  Cesium 

I 

Standard  Deviation  -  Element  C  above 

J 

"  "  -  Element  D  above 

K 

L 

”  "  -  Element  E  above 

"  "  -  Element  F  above 

-  Element  G  above 

-  Element  H  above 

GPS  -  Online  Cesium  Measurement  Residual 
GPS  -  Secondary  Cesium  Measurement  Residual 
GPS  -  Primary  Cesium  Measurement  Residual 
Online  -  Secondary  Cesium  Measurement  Residual 
Online  -  Primary  Cesium  Measurement  Residual 


M 

N 

O 

P 

Q 

R 

S 


CONTENT 


ELEMENT 
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0.05323 

93.67398 

43.43991  -27.6859  -0.67432 
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13 
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25.19416 
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13 

47605 

9  -8.72544 

0.00974 

32.82698 

-0.3704  -27.6606  -0.27409 
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0.00974 

32.64178 
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Figure  4.4-1  Sample  Worksheet  Generated  from  LOTUS.OUT  File 


MENU 


Worksheet  Range  Copy 

Move  File  Print 

Graph 

Data  System  Add-In  Quit 

Global  Insert 

Delete 

Column 

Erase  Titles  Window  Status  Page  Learn 

I 

J 

K 

L 

M 

N 

O  P 

1 

0.00981 

0.04 

0.01437 

0.04 

0 

0.04 

-0.26809  -104.602 

0. 02926 

0.04037 

0.03109 

0.06782 

0.02757 

0.04037 

-0.26809  -104.602 

0 . 04923 

0.04074 

0.05724 

0.08718 

0.04824 

0.04074 

-0.26809  -104.602 

3 .06938 

0.04111 

0.09293 

0.10296 

0.06868 

0.04111 

-0.26809  -104.602 

z, 

0.08967 

0.04147 

0.13619 

0.11662 

0.08913 

0.04147 

-0.26809  -104.602 

6 

0 . 1101 

0.04183 

0.18576 

0.12884 

0.10966 

0.04183 

-0.26809  -104.602 

0.0467 

0.02213 

0.04898 

0.08535 

0.04669 

0.02198 

-0.16761  -96.0366 

3 

0.05749 

0.0228 

0.07555 

0.10141 

0.0575 

0.02266 

-0.16761  -96.0366 

9 

0.06843 

0.02345 

0.1158 

0.11526 

0.06842 

0.02331 

-0.16761  -96.0366 

10 

0 . 07954 

0.02408 

0.16433 

0.12761 

0 . 0795 

0.02394 

-0.16761  -96.0366 

11 

0. 09085 

0.02469 

0.21907 

0.13887 

0.09077 

0.02456 

-0.16761  -96.0366 

12 

0.10236 

0.02529 

0.27906 

0.14928 

0.10224 

0.02517 

-0.16761  -96.0366 

13 

0 . 0458 

0.01783 

0.0495 

0.09021 

0.04579 

0.01776 

0.13827  126.0512 

14 

0 . 05398 

0.01865 

0 . 07608 

0.10554 

0.05395 

0.01859 

0.13827  126.0512 

15 

0.06242 

0.01944 

0.11809 

0.1189 

0.06237 

0.01938 

0.13827  126.0512 

16 

0.07113 

0.0202 

0.16858 

0.13091 

0.07105 

0.02013 

3.13827  126.0512 

0  .  08013 

3.02093 

3.2252 

0.14191 

0 . 08002 

3.02087 

J.1JS2/  ...  6  .  05  I* 

0  .  C894 

0.02163 

0.28693 

0.15211 

0.08926 

3.02157 

0.13327  1 2  6  .  w  5 1 .. 

19 

0 . 04463 

0.01735 

0.04953 

0.09023 

0.04462 

3 . 01733 

0.17911  -9.86731 

20 

I  .  C  5  2  3  8 

0.0182 

0.07609 

0. 1C555 

0.05236 

3. 01818 

3.17911  -9. 8673. 

Figure  4.4-1  (con’t)  Elements  I  through  P  of  the  Sample  Worksheet 
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Figure  4.4-1  (con’t)  Elements  Q  through  S  of  the  Sample  Worksheet 
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TEST  RESULTS 


This  chapter  presents  the  results  of  several  scenarios  examined  using  the  CSP.  The  main 
purpose  of  these  runs  is  to  examine  the  steady  state  behavior  of  the  filter  in  different  configurations. 
The  section  focuses  on  the  data  collected  at  the  OMSTA  Norway.  The  RDAS  timing  data  used 
in  this  effort  was  provided  to  SYNETICS  by  USNO  personnel;  specific  data  intervals  provided 
were: 

•  OMSTA  Norway:  11  March  1989  (47596)  through  12  April  1989  (47628) 

•  OMSTA  Hawaii:  8  March  1989  (47593)  through  12  April  1989  (47628) 

The  results  presented  in  Section  5.1  were  developed  using  a  subset  of  this  data.  Additional  plots 
for  both  the  Norway  and  Hawaii  Omega  Stations  are  provided  in  a  separate  addendum  (Reference 
2). 

Before  presenting  the  data  plots,  a  brief  review  of  the  thought  process  underlying  this  study 
will  be  useful.  As  mentioned  in  Chapter  3,  the  data  collected  by  USNO  is  independent  of  the 
current  Omega  synchronization  process.  Specifically,  the  cesium  timing  is  measured  before  any 
|  OMSFOG  corrections.  This  allows  us  to  use  the  collected  data,  without  having  to  "back  out" 
corrections,  to  hypothesize  a  new  control  scheme.  An  important  part  of  the  current  study  is  to 
investigate  how  frequently  corrections  should  be  made  in  any  new  control  scheme. 

To  proceed  with  this  investigation,  we  will  simulate  obtaining  data  twice  a  day  and 
estimating/correcting  at  the  following  rates:  once  a  day,  once  every  third  day,  once  every  fifth  day, 
and  once  every  seventh  day.  A  comparison  of  the  resulting  performances  will  give  an  indication 
of  benefits  (which  must  be  considered  within  the  context  of  the  entire  Omega  system  error  budget) 
to  be  derived  from  more  frequent  corrections. 

For  any  given  configuration,  a  correction  or  measurement  update  interval  is  selected. 
Although  we  say  data  is  "collected"  twice  a  day,  we  simulate  an  implementation  scheme  in  which 
it  is  only  provided  as  an  input  to  the  estimation  filter  at  prescribed  intervals  (daily,  every  third  day, 
etc.)  The  output  of  each  measurement  update  filter  cycle  (performed  at  the  prescribed  interval) 
is  an  estimate  of  the  phase  and  frequency  offsets  of  the  cesiums  and  the  corresponding  filter- 
indicated  RMS  errors,  all  valid  at  the  measurement  update  time.  To  obtain  estimates  at  times  in 
between  measurement  incorporation  events,  the  phase  offset  estimates  are  extrapolated  forward 


to  the  desired  time  using  the  frequency  offset  estimate.  When  the  end  of  the  prescribed 
measurement  update  interval  is  reached,  this  extrapolated  estimate  is  corrected  with  the  new 
measurement  data. 

The  degree  to  which  the  extrapolated  phase  offset  estimate  accurately  reflects  the  true 
phase  offset  will  clearly  decrease  as  the  measurement  update  interval  is  extended.  A  measure  of 
this  deviation  is  provided  by  the  phase  offset  residual  plots.  These  plots  show  the  difference 
between  the  predicted  phase  offsets  at  the  time  at  which  a  set  of  corrections  is  requested,  and  the 
actual  measurement  data  which  is  "revealed  to  the  filter"  an  "update  interval"  later.  An  assessment 
of  these  residuals,  for  various  update  intervals,  is  what  constitutes  the  major  output  of  this  study. 

5.1  TEST  CONFIGURATIONS 

Five  distinct  measurement  sets  are  examined  with  each  set  defining  a  configuration.  Within 
each  configuration,  the  update  rates  were  varied  with  a  maximum  update  interval  of  7  days.  The 
selected  configurations  are  illustrated  pictorially  in  Figure  5.1*1.  For  each  configuration,  Figure 
5.1-1  shows  the  four  main  units:  the  online  cesium,  the  secondary  cesium,  the  primary  cesium,  and 
the  GPS  Datum  Monitor,  with  directed  arrows  joining  them.  The  arrows  designate  which 
measurements  are  defined  in  the  configuration.  Configuration  I  is  discussed  in  some  detail  with 
the  intent  of  explaining  the  significance  of  the  various  plot  types  provided  by  the  CSP.  The 
measurement  sets  available  at  the  OMSTAs  Hawaii  and  Norway  are  provided  in  Table  5.1-1. 
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5.1-1  Test  Configurations 


TABLE  5.1-1 

RDAS  TIMING  DATA  MEASUREMENTS 


PHASE  OFFSET 

MEASUREMENT 

OMEGA 

STATION 

MEASUREMENT  DESIGNATION 

Time  difference  between  GPS 

and  online  cesium 

Norway 

GPS-CS486=(GPS-CS124)-(CS486-CS124) 

Hawaii 

GPS-CS529=(GPS-CS554)-(CS529-CS554) 

Time  difference  between  GPS 

and  secondary  cesium 

Norway 

GPS-CS017 

Hawaii 

GPS-CS349 

Time  difference  between  GPS 

and  primary  cesium 

Norway 

GPS-CS124 

Hawaii 

GPS-CS554 

Time  difference  between 

online  and  secondary  cesiums 

Norway 

CS486-CS017 

Hawaii 

CS529-CS349 

Time  difference  between 

online  and  primary  cesiums 

Norway 

CS486-CS124 

Hawaii 

CS529-CS554 
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5.1.1  Configuration  I  (Norway) 


In  configuration  I  the  available  measurement  set  was  assumed  to  consist  of  GPS/cesium 
measurements  for  each  of  the  three  cesiums:  online,  primary,  and  secondary.  With  this 
measurement  set,  the  six-state  filter  reduces  to  a  set  of  three  uncoupled  two  state  filters.  In  this 
uncoupled  configuration  the  estimate  of  the  online  clock  phase  offset  and  frequency  depend  only 
on  the  GPS/  online  cesium  measurements.  This  same  relation  between  the  measurements  and  the 
state  estimates  also  exists  for  the  primary  and  secondary  clocks. 

Several  different  runs  were  performed  in  this  configuration,  with  the  interval  between 
measurement  updates,  that  is  the  interval  at  which  measurements  are  provided  to  the  filter,  varied 
between  0.5  days  and  7  days.  In  the  following  two  subsections  the  0.5  day  and  5  day  interval  cases 
are  examined.  Since  the  run  with  the  0.5  day  update  interval  is,  within  the  context  of  the  cases 
examined,  a  best  case,  this  run  will  serve  as  a  baseline  with  which  to  compare  the  results  of  other 
runs  in  this  configuration  as  well  as  runs  in  other  configurations. 

5. 1.1.1  Configuration  /  -  0.5  day  update  interval 

The  filter’s  phase  and  frequency  estimates  for  each  of  the  three  cesiums  are  provided  in 
Figures  5.1-2  and  5.1-3,  respectively.  Note  that  for  plotting  purposes,  the  initial  phase  offset 
estimate  used  by  the  filter  is  subtracted  cut  of  the  data  record  ensuring  each  of  the  records  begin 
at  zero.  The  online  clock  (CS486)  phase  offset  estimate  is  seen  to  remain  very  stable  throughout 
the  interval.  This  is  consistent  with  the  frequency  offset  estimate  for  CS486  which  remains  near 
zero.  Both  CS124  and  CS017  show  more  evident  phase  offset  drifts.  The  frequency  offset 
estimates  for  these  two  cesiums  are  somewhat  less  stable  than  those  of  CS486. 

The  phase  offset  residuals,  given  in  Figure  5.1-4,  provide  a  measure  of  the  performance  of 
the  filter.  For  this  minimum  update  rate  case,  these  residuals  reduce  to  the  difference  between  the 
actual  measurement  and  the  measurement  updated  estimate.  Throughout  the  interval  these 
residuals  remain  below  the  0.1  microsecond  (100  nanosecond)  level.  The  largest  residuals  are  those 
of  CS124  which  occur  near  day  20.  As  is  evident,  a  high  degree  of  correlation  exists  between  the 
three  phase  offset  residuals.  The  noise  in  each  of  the  three  GPS/cesium  phase  measurements  is 
mainly  the  result  noise  in  the  GPS  outputs,  and  is  thus  common  to  each  of  the  three  measurements. 
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The  filter-indicated  RMS  errors  for  each  of  the  three  phase  offset  states  are  provided  in 
Figure  5.1-5.  These  values  indicate  the  filter’s  estimate  of  the  uncertainty  in  the  provided  estimates. 
Although  the  true  errors  are  not  available  for  comparison,  the  relative  consistency  of  the  phase 
offset  residuals  with  the  filter  indicated  RMS  values  is  an  indication  of  reasonable  performance. 
A  large  discrepancy  between  these  two  quantities,  specifically  a  diverging  residual  and  stable  filter 
RMS  value,  may  be  taken  as  an  indication  of  mismodeling.  No  problem  of  this  type  is  evident 
here. 


5. 1.1.2  Configuration  /  -  5  day  update  interval 

In  this  run  the  same  configuration  as  above  was  used  but  the  interval  between  measurement 
updates  was  lengthened  from  0.5  days  to  5  days.  Plots  summarizing  this  run  are  provided  in  Figures 
5.1-6  through  5.1-9.  The  updated  phase  estimates  from  day  one  of  an  update  interval  are 
propagated  forward  using  the  updated  frequency  estimates  until  the  next  measurement  becomes 
available  to  the  filter.  When  this  next  measurement  is  received,  the  filter’s  measurement  update 
step  is  executed  and  the  propagated  estimate  is  updated.  This  propagation  and  update  cycle  is 
evident  in  the  state  estimates  (Figure  5.1-6)  and  the  filter  RMS  error  estimate  (Figure  5.1-9).  The 
measurement  update  events  are  evident  in  these  plots  occurring  at  days  5,  10,  etc. 

The  phase  offset  estimates  show  the  same  general  behavior  as  was  observed  in  the  0.5  day 
update  interval  case.  The  frequency  offset  estimates  are  also  similar  to  those  of  the  0.5  day  update 
interval  case.  Since  no  frequency  rate  terms  are  maintained  by  the  filter,  the  frequency  offset 
estimates  are  held  constant  between  measurement  updates  and  change  only  at  these  events. 

The  phase  offset  residuals  (Figure  5.1-8)  again  provide  a  measure  of  the  error  in  the  phase 
offset  estimates.  The  residuals  for  all  three  clocks  remain  below  0.25  microseconds  (250 
nanoseconds)  throughout  the  interval.  The  largest  residuals  are  those  associated  with  CS124  and 
occur  near  day  24  and  day  29,  both  days  corresponding  to  ends  of  a  filter  update  cycle  with 
measurement  updates  occurring  on  days  25  and  30.  Some  portion  of  this  error  may  be  attributed 
to  an  apparent  error  in  the  frequency  offset  estimates  for  CS124  near  these  days.  An  error  in  the 
frequency  offset  estimate  on  the  order  of  0.03  microseconds/day,  for  example,  will  propagate  into 
a  phase  offset  estimate  error  of  magnitude  0.15  microseconds  at  the  end  of  a  five  day  period.  Some 
measure  of  the  eiror  in  the  frequency  offset  estimate  in  this  case  may  be  gained  by  comparing  with 
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the  corresponding  frequency  offset  for  the  0.5  day  update  interval  case.  At  the  measurement 
update  times  (days  5,  10,  15,  etc.)  the  phase  offset  residuals  for  each  of  the  clocks  drops  to  near 
zero.  The  filter  has  effectively  "reset"  its  phase  offset  estimate  to  the  value  of  the  newly  received 
phase  offset  measurement.  This  behavior  is  consistent  with  the  relative  smallness  of  the  assumed 
uncertainty  in  the  GPS/cesium  measurements  as  compared  to  the  uncertainty  in  the  a  priori 
estimate,  as  expressed  in  the  propagated  error  covariance. 

Much  of  the  noise  evident  in  the  each  of  the  three  residuals  is  a  result  of  noise  in  the  GPS/ 
cesium  measurements,  rather  than  in  the  estimates  themselves.  This  may  be  inferred  since  the 
estimates  change  smoothly  (linearly)  between  updates,  and  may  also  be  verified  by  examining  the 
raw  data.  As  in  the  0.5  day  case,  a  high  degree  correlation  exists  between  the  three  residuals. 

The  filter  indicated  RMS  errors  (Figure  5.1*9)  range  between  0.05  microseconds,  just  after 
a  measurement  update  and  0.13  microseconds  just  before  an  update.  The  growth  between 
estimates  is  due  to  the  inclusion  of  the  process  noise  matrix.  Note  that  the  post-update  RMS  error 
(0.05  microseconds)  matches  the  value  assumed  for  the  GPS/cesium  measurements  (aGC  =  0.05 
|  fisec). 


5.1.2  Configuration  II  (Norway) 

As  in  configuration  I,  the  measurement  set  assumed  here  is  sufficient  to  obtain  reasonable 
estimates  of  the  phase  offsets  for  each  of  the  cesiums.  In  this  case,  however,  only  one  direct 
GPS/cesium  measurement  is  used.  To  complete  the  measurement  set,  online/primary  and 
online/secondary  measurements  are  included.  The  0.5  day  update  interval  case  is  summarized  in 
Figures  5.1-10  through  5.1-13  and  the  5  day  update  interval  case  is  summarized  in  Figures  5.1-14 
through  5.1-17.  The  basic  behavior  of  the  phase  and  frequency  offset  plots  is  the  same  as  was 
observed  for  Configuration  I.  Some  slight  differences  in  the  residual  plots  for  the  two 
configurations  can  be  seen.  The  phase  offset  residuals  for  Configuration  II  tend  to  be  slightly  larger 
when  compared  to  the  same  update  interval  case  in  Configuration  I.  Larger  phase  offset  residuals 
for  the  primary  and  secondary  cesiums  might  be  expected,  since  only  indirect  measurement  of  the 
phase  offsets  of  these  clocks  are  used  in  this  configuration. 
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5.2 


DISCUSSION 


In  this  chapter  CSP  results  for  configurations  I  and  II  have  been  presented  using  Norway 
data.  These  two  configurations  represent  nominal  cases  in  that  a  sufficient  measurement  set  is 
available  to  estimate  the  states  of  each  of  the  three  clocks.  For  each  of  the  two  configurations,  two 
different  measurement  update  schedules,  0.5  days  and  5  days,  were  examined.  These  runs  served 
to  confirm  the  integrity  of  the  basic  filter  design.  In  all  cases  reasonable  estimates  were  maintained 
and  no  instances  of  divergent  behavior  were  noted.  No  significant  difference  between  the  two  main 
configurations  were  observed,  although  the  phase  offset  residuals  for  the  configuration  I  runs  were 
slightly  smaller  in  magnitude.  The  largest  phase  offset  residual  recorded  during  the  runs  presented 
in  this  chapter  was  on  the  order  of  300  nanoseconds. 

The  addendum  that  accompanies  this  report  (Reference  2)  examines  additional 
measurement  configurations  and  measurement  update  intervals,  using  both  Hawaii  and  Norway 
data.  In  reviewing  these  plots,  the  reader  should  focus  his  attention  on  the  phase  offset  residual 
plots.  To  help  in  a  review  of  this  data,  a  record  of  the  maximum  residual  recorded  for  each  of  the 
three  configurations  involving  GPS/cesium  measurement  (configurations  I,  II,  and  III)  was  recorded. 
This  information  was  obtained  by  reviewing  the  phase  offset  residual  plots  for  each  run  and 
selecting  the  absolute  maximum  value  and  is  provided  in  Table  5.2-1  for  OMSTA  Norway  and 
Table  5.2-2  for  OMSTA  Hawaii.  In  configuration  III  the  measurement  set  is  such  that  the  phase 
offset  estimate  cannot  be  computed  for  one  of  the  cesiums  (i.e.,  CS124  at  Norway,  and  CS554  at 
Hawaii);  the  corresponding  columns  in  the  tables  are  thus  left  blank.  Note  the  variations  between 
the  different  configurations  and  different  update  intervals  is  relatively  small.  This  is  primarily  due 
to  the  fact  that  the  frequency  offset  estimates  of  the  cesiums  remain  relatively  constant  over  the 
interval. 

The  same  approach  was  taken  with  the  Hawaii  data  as  was  used  with  the  Norway  with  one 
exception.  Due  to  GPS  coverage  limitations  at  the  Hawaii  station.  GPS  measurement  data  was  not 
always  available  at  each  half  day  interval  as  was  the  case  at  Norway.  The  topography  at  OMSTA 
Hawaii  is  such  that  satellites  at  low  elevations  are  obscured  for  a  large  range  of  azimuths.  This 
fact,  combined  with  the  fact  that  only  a  limited  constellation  was  in  place,  resulted  in  the  periodic 
loss  of  GPS  measurement  data.  To  overcome  this  limitation  in  the  CSP  runs,  the  basic  minimum 
interval  at  which  measurements  were  used  for  residual  computations  in  the  CSP  was  increased  to 
1  day  for  the  Hawaii  scenarios.  This  spacing  was  found  to  be  sufficient  to  avoid  the  measurement 
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outages.  In  the  Norway  scenarios  examined  in  this  report,  and  provided  in  the  addendum,  the  basic 
interval  of  1/2  day  was  used.  To  help  the  reader  interpret  the  Hawaii  plots  provided  in  the 
addendum,  a  sample  Hawaii  phase  offset  residual  plot  is  provided  in  Figure  5.2-1  corresponding 
to  the  five  day  measurement  update  interval  scenario.  This  plot  shows  the  residual  computed  at 
1  day  intervals,  rather  than  1/2  day  intervals. 

The  initial  phase  offset  estimates  for  the  cases  examined  above,  as  well  as  for  those 
examined  in  the  addendum,  were  set  at  the  initial  GPS/cesium  measurements.  This  selection  was 
made  since  our  main  interest  is  in  the  steady  state  response  of  the  filter.  This  aspect  of 
performance  is  the  most  relevant  for  any  hypothesized  synchronization  scheme  that  uses  GPS 
data.  Although  of  less  interest  for  the  objectives  of  this  report,  the  transient  response  of  the  filter 
was  also  examined  for  additional  verification  of  filter  performance.  This  was  accomplished  by 
injecting  an  arbitrary  error  in  the  initial  phase  offset  estimate  and  verifying  that  the  estimates 
eventually  converged  to  reasonable  values.  A  sample  case  that  illustrates  this  transient  behavior 
is  provided  in  Figures  5.2-2  and  5.2-3,  showing  the  phase  offset  estimates  and  phase  residuals, 
respectively.  The  initial  errors  injected  for  this  test  case  ranged  between  0.3  and  1.2  microseconds. 

)  Although  a  longer  time  period  is  required  for  the  filter  to  reach  the  point  were  it  is  maintaining 
stable  estimates,  this  point  is  eventually  reached. 
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TABLE  5.2-1 

NORWAY  OMEGA  STATION  RESIDUALS 


TABLE  5.2-2 

HAWAII  OMEGA  STATION  RESIDUALS 


ABSOLUTE  MAXIMUM  PHASE  OFFSET  RESIDUALS  (microseconds) 


Configuration  I 

CS529 

CS349 

CS554 

1.0  day  update 

.053 

.084 

.092 

3.0  day  update 

.115 

.159 

.185 

5.0  day  update 

.148 

.189 

.154 

7.0  day  update 

.178 

.189 

.120 

Configuration  11 

CS529 

CS349 

CS554 

1.0  day  update 

.069 

.109 

.101 

3.0  day  update 

.112 

.158 

.173 

5.0  day  update 

.147 

.189 

.149 

7.0  day  update 

.154 

.189 

.134 

Configuration  III 

CS529 

CS349 

CS554 

1.0  day  update 

.058 

.099 

3.0  day  update 

.109 

.152 

5.0  day  update 

.148 

.189 

7.0  day  update 

.170 

.189 
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Figure  5.1-2  Configuration  I^ftase  Offset  Estimates  (0.5  day) 
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Figure  5.1-3  Configuration  I  -  Frequency  Offset  Estimates  (0.5  day) 
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Figure  5.1-4  Configuration  l  -  Phase  Offset  Residuals  (0.5  day) 
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Figure  5.1-6  Configuration  I  -  Phase  Offset  Estimates  (5  day) 
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Figure  5.1-7  Configuration  I  -  Frequency  Offset  Estimates  (5  day) 
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Figure  5.1-8  Configuration  I  -  Phase  Offset  Residuals  (5  day) 
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F  igure  5.1-9  Conliguration  I  -  Filler  Indicated  RMS  Phase  Offset  Errors  (5  day) 
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Figure  5.1-1 1  Configuration  II  -  Frequency  Offset  Estimates  (0.5  day) 
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12  Configuration  II  -  Ph^ Offset  Residual  (0.5  day) 
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Hgure  5.1-13  Configuration  II  -  Filter  Indicated  RMS  Phase  Offset  Errors  (0.5  day) 
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5.1-14  Configuration  II  -  Offset  Estimates  (5  day) 
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5.1-15  Configuration  H  -  Frequency  Offset  Estimates  (5  day) 
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1. 


INTRODUCTION 


To  demonstrate  the  Covariance  Simulation  Program  capabilities,  five  different 
configurations  were  used  with  data  from  Omega  Stations  Norway  and  Hawaii  Chapters  2  and 
3  contain  the  resulting  plots.  The  organization  and  content  of  these  chapters  are  described  below. 

Chapter  2  begins  with  the  CSP  test  configuration  for  Norway  Station  followed  by  five 
sections  containing  the  graphs  for  the  corresponding  configuration.  Each  section  is  identified  by 
the  configuration  number  and  includes  an  alphabetical  list  of  the  plots  generated  for  the  particular 
configuration.  All  of  the  configurations  were  executed  at  1  day,  3  day,  5  day  and  7  day  update 
rates.  Chapter  3  has  an  identical  layout  but  was  generated  on  a  month  of  data  for  the  Hawaii 
Station.  Chapter  4  includes  a  source  code  listing  of  the  main  Kalman  filter  routine,  GPSSIM,  and 
the  file  translation  routine,  FILEPRC. 


2. 


NORWAY 


Table  2-1  depicts  five  CSP  configurations  that  were  used  to  process  RDAS  timing  data 
from  Omega  Station  Norway. 


Thble  2-1 

CSP  TEST  CONFIGURATION 
(Norway  Station) 


CONFIGURATION 

NO. 

GPS  TIME 
TRANSFER 

ON-LINE 

CESIUM 

(CS486) 

PRIMARY 

CESIUM 

(CS124) 

SECONDARY 

CESIUM 

(CS017) 

I 

On 

On 

On 

On 

II 

On 

On 

On 

On 

III 

On 

On 

Off 

On 

IV 

Off 

On 

On 

On 

V 

Off 

On 

Off 

On 

2- 


2.1 


CONFIGURATION  I 


Seven  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order: 

a)  Phase  Offset  Estimate 

b)  Frequency  Offset  Estimate 

c)  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017,  GPS-CS124 

d)  Phase  Covariance 

e)  Frequency  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPS1,  GPS2 
and  GPS3.  Event  "GPS1"  identifies  GPS-CS486  measurement;  event  "GPS2"  identifies  GPS-CS017 
measurement  and  event  "GPS3"  identifies  GPS-CS124  measurement.  The  remainder  of  this 
section  contains  plots  as  described  above. 
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Figure  2.1-lb  Frequency  Offset  Estimate  (1  Day  Update) 
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RMS  ESTIMATION  ERROR  Imico 
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Figure  2.1-lc  Phase  Offset  Residual:  GPS-CS4«6,-CS017,-CS124  (1  Day  Update) 
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Figure  2.1-le  Frequency  Covariance  (1  Day  Update) 
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Figure  2.1-2e  Frequency  Covariance  (3  Day  Update) 
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GPS/ OMEGA  COVARIANCE  SIMULATION 
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Figure  2.1-3e  Frequency  Covariance  (5  Day  Update) 
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Figure  2.1-4c  Phase  Offset  Residual:  GPS-CS486,-CS017,-CS124  (7  Day  Update) 
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GPS/OMEGA  COVARIANCE  SIMULATION 
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Figure  2.1-4e  Frequency  Covariance  (7  Day  Update) 
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2.2 


CONFIGURATION  U 


Five  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the  following 
order: 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017,  GPS-CS124 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPSl,  CLKI 
and  CLK2.  Event  "GPSl”  identifies  GPS-CS486  measurement;  event  "CLKI"  identifies  CS486- 
CS017  measurement  and  event  "CLK2"  identifies  CS486-CS124  measurement.  The  remainder  of 
this  section  contains  plots  as  described  above. 
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GPS/OMEGA  COVARIANCE  SIMULATION 
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Figure  2.2<3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  2.2-3b  Phase  Offset  Residual:  GPS-CS486,.CS017,-CS124  (5  Day  Update) 
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Figure  22-4*  Phase  Offset  Estimate  (7  Day  Update) 
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Figure  22-4b  Phase  Offset  Residual:  GPS-CS486,*CS017  (7  Day  Update) 
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Figure  12-4e  Phue  Covarteace  (7  Day  Update) 
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CONFIGURATION  m 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPSl  and  CLK1. 
Event  "GPSl"  identifies  GPS-CS486  measurement  and  event  "CLK1"  identifies  CS486-CS017 
measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  23-la  Phase  Offset  Estimate  (1  Day  Update) 
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Figure  23-lb  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (1  Day  Update) 
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Figure  23-lc  Phase  Covariance  (1  Day  Update) 
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Figure  23-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  2-3-3b  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (5  Day  Update) 
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Figure  23-3c  Phase  Covariance  (5  Day  Update) 
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Figure  23-4b  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (7  Day  Update) 
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Figure  2-3-4c  Phase  Covariance  (7  Day  Update) 
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CONFIGURATION  IV 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order: 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017,  GPS-CS124 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  CLK1  and  CLK2. 
Event  "CLK1"  identifies  CS486-CS017  measurement  and  event  "CLK2"  identifies  CS486-CS124 
measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  2.4-lc  Phase  Covariance  (1  Day  Update) 
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Figure  2.4-2a  Phase  Offset  Estimate  (3  Day  Update) 
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Figure  2.4-2b  Phase  Offset  Residual:  GPS-CS486,  -CS017,  -CS124  (3  Day  Update) 
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Figure  2.4-2c  Phase  Covariance  (3  Day  Update) 
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Figure  2.4-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  2.4-3b  Phase  Offset  Residual:  GPS-CS486,  -CS017,  -CS124  (5  Day  Update) 
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Figure  2.4-3c  Phase  Covariance  (5  Day  Update) 
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Figure  2.4-4a  Phase  Offset  Estimate  (7  Day  Update) 
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CONFIGURATION  V 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order: 

a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017 

c)  Phase  Covariance 

This  configuration  was  implemented  by  the  following  event:  CLK1.  This  event  identifies 
CS486-CS017  measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  2.5-lb  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (1  Day  Update) 
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Figure  2.5-lc  Phase  Covariance  (1  Day  Update) 
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Figure  2.5-2a  Phase  Offset  Estimate  (3  Day  Update) 
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Figure  2J«2b  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (3  Day  Update) 
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Figure  2.5-2c  Phase  Covariance  (3  Day  Update) 
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Figure  2.5-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  2.5*3b  Phase  Offset  Residual:  GPS-CS486,  GPS-CS017  (5  Day  Update) 
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Figure  2.5-3c  Phase  Covariance  (5  Day  Update) 
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Figure  2-5-4c  Phase  Covariance  (7  Day  Update) 
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3. 


HAWAII 


Table  3-1  depicts  five  CSP  configurations  that  were  used  to  process  RDAS  timing  data  from 
Omega  Station. 


Thble  3-1 

CSP  TEST  CONFIGURATION 
(Hawaii  Station) 


CONFIGURATION 

NO. 

GPS  TIME 
TRANSFER 

ON-LINE 

CESIUM 

(CS529) 

PRIMARY 

CESIUM 

(CS554) 

SECONDARY 

CESIUM 

(CS349) 

I 

On 

On 

On 

On 

II 

On 

On 

On 

On 

III 

On 

On 

On 

Off 

IV 

Off 

On 

On 

On 

V 

Off 

On 

On 

Off 

3.1 


CONFIGURATION  I 


Seven  types  of  plots  were  generted  for  this  configuration.  They  are  presented  in  the  following 
order: 


a)  Phase  Offset  Estimate 

b)  Frequency  Offset  Estimate 

c)  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349 

d)  Phase  Covariance 

e)  Frequency  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPS1,  GPS2 
and  GPS3.  Event  "GPS1"  identifies  GPS-CS529  measurement;  event  "GPS2"  identifies  GPS-CS554 
measurement  and  event  "GPS3"  identifies  GPS-CS349  measurement.  The  remainder  of  this 
section  contains  plots  as  described  above. 
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Figure  3.1-la  Phase  Offset  Estimate  (1  Day  Update) 

GPS/OMEGA  COVARIANCE  SIMULATION 

FRED  OFFSET  ESTMATE  (HUT  UPDATE) 


0  I  2  :  4  I  6  I  •  l  10  I  12  I  14  I  W  i  18  I  20  I  22  1  24  i  26  I 

1  3  5  7  9  11  13  IS  17  19  21  23  25  2i 

a  CSS29  -  CS5S4  o  CS34S 

Figure  3.1-lb  Frequency  Offset  Estimate  (1  Day  Update) 
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Figure  3.1-lc  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (1  Day  Update) 
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Figure  3.1-ld  Phase  Covariance  (1  Day  Update) 
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GPS/ OMEGA  COVARIANCE  SIMULATION 


Figure  3.1-2c  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (3  Day  Update) 
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Figure  3.1*2d  Phase  Covariance  (3  Day  Update) 
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Figure  3.1-2e  Frequency  Covariance  (3  Day  Update) 
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Figure  3.1*3b  Frequency  Offset  Estimate  (5  Day  Update) 
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Figure  3.1-3c  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (5  Day  Update) 
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Figure  3.1-3d  Phase  Covariance  (5  Day  Update) 
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Figure  3.1-3e  Frequency  Covariance  (5  Day  Update) 
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Figure  3.1-4c  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (7  Day  Update) 
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CONFIGURATION  II 


Five  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the  following 
order: 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPS1,  CLK1 
and  CLK2.  Event  "GPSl"  identifies  GPS-CS529  measurement;  event  "CLKl"  identifies  CS529- 
CS554  measurement  and  event  "CLK2"  identifies  CS529-CS349  measurement.  The  remainder 
of  this  section  contains  plots  as  described  above. 
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Figure  3.2-lb  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (1  Dav  Update) 
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Figure  3.2-lc  Phase  Covariance  (1  Day  Update) 
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Figure  3.2-2a  Phase  Offset  Estimate  (3  Day  Update) 
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Figure  3.2-2b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (3  Day  Update) 
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Figure  32-2c  Phase  Covariance  (3  Day  Update) 
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Figure  3.2-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  3.2-3b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (5  Day  Update) 
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Figure  3.2-3c  Phase  Covariance  (5  Day  Update) 
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Figure  3.2-4a  Phase  Offset  Estimate  (7  Day  Update) 
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Figure  3.2-4b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (7  Day  Update) 
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Figure  3.2-4c  Phase  Covariance  (7  Day  Update) 
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3.3 


CONFIGURATION  HI 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order. 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  types  of  events:  GPS1  and  CLK1. 
Event  "GPS1"  identifies  GPS-CS529  measurement  and  event  "CLK1"  identifies  CS529-CS554 
measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  3J-2a  Phase  Offset  Estimate  (3  Day  Update) 
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Figure  3.3-2b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554  (3  Day  Update) 
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Figure  3  J-2c  Phase  Covariance  (3  Day  Update) 
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Figure  33-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  33-3b  Phase  Offset  Residual:  GPS-CS529,  GPS-CSS54  (5  Day  Update) 
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Figure  33-3c  Phase  Covariance  (5  Day  Update) 
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Figure  3 .3 -4a  Phase  Offset  Estimate  (7  Day  Update) 
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Figure  33-4b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554  (7  Day  Update) 
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Figure  3 .3 -4c  Phase  Covariance  (7  Day  Update) 
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3.4 


CONFIGURATION  IV 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order: 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349 

c)  Phase  Covariance 


i 


This  configuration  was  implemented  by  the  following  types  of  events:  CLKl  and  CLK2. 
Event  "CLKl"  identifies  CS529-CS554  measurement  and  event  "CLK2"  identifies  CS529-CS349 
measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  3.4-la  Phase  Offset  Estimate  (1  Day  Update) 
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Figure  3.4-lb  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (1  Day  Update) 
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Figure  3.4- lc  Phase  Covariance  (1  Day  Update) 
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Figure  3.4-2b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (3  Day  Update) 
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Figure  3.4-2c  Phase  Covariance  (3  Day  Update) 
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Figure  3.4-3a  Phase  Offset  Estimate  (5  Day  Update) 
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Figure  3.4-3b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554,  GPS-CS349  (5  Day  Update) 
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Figure  3.4-3c  Phase  Covariance  (5  Day  Update) 
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Figure  3.4-4a  Phase  Offset  Estimate  (7  Day  Update) 
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Figure  3.4-4b  Phase  Offset  Residual:  GPS-CSS29,  GPS-CS554,  GPS-CS349  (7  Day  Update) 


3-40 


0  1  ^  1  *  i  ®  1  ®  I  10  I  1J  l  14  I  18  l  I#  i  20  |  22  I  24  I  26  . 

*  3  s  r  9  H  13  IS  17  19  21  23  25  27 

TMBdayai 

□  CS529  *  CSS54  o  CS349 


Figure  3.4-4c  Phase  Covariance  (7  Day  Update) 
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3.5 


CONFIGURATION  V 


Three  types  of  plots  were  generated  for  this  configuration.  They  are  presented  in  the 
following  order: 


a)  Phase  Offset  Estimate 

b)  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554 

c)  Phase  Covariance 


This  configuration  was  implemented  by  the  following  event:  CLK1.  This  event  identifies 
CS529-CS554  measurement.  The  remainder  of  this  section  contains  plots  as  described  above. 
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Figure  3.5-lb  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554  (1  Day  Update) 
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Figure  3.5-lc  Phase  Covariance  (1  Day  Update) 
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Figure  3.5-2a  Phase  Offset  Estimate  (3  Day  Update) 
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Figure  3.5-2b  Phase  Offset  Residual:  GPS-CS529,  GPS-CS554  (3  Day  Update) 
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Figure  3.5-2c  Phase  Covariance  (3  Day  Update) 
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Figure  3.5-3a  Phase  Offset  Estimate  (5  Day  Update) 

GPS/OMEGA  COVARIANCE  SIMULATION 


PHASE  OFFSET  RESBUAL  15-OAT  IPOATF1 


THtCldayai 

a  GPS-CSS29  -  CPS-CSSM 


Figure  3.5-3b  Phase  Offset  Residual:  GPS-CS529,  GPS-CSS54  (5  Day  Update) 
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Figure  3.5-3c  Phase  Covariance  (5  Day  Update) 
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Figure  3.5-4c  Phase  Covariance  (7  Day  Update) 
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C 

C 

C 

C 


C 

C 

C 

C 

C 

C 

C 


101 


1 


102 


10 


1 


105 


PROGRAM  FILEPRC 

THIS  PROGRAM  CONVERTS  RAW  INPUT  DATA  FILE  SPECIFIED  BY 

INFILE  INTO  SUITABLE  FORMAT  ( NORWAY . DAT )  FOR  GPSSIM  (FILTER)  PROGRAM 

CHARACTER* 1  INDAT(80) 

CHARACTER* 15  INFILE, OUTFILE 

R£AL*8  DUMMY, GPSCS2  ,GPSCS3  ,  CS1CS2  ,  CS1CS3 

REAL*  8  MJDAY , THRESH 

INTEGER  SYSTIME 


GPSCS2  —>  GPS  -  CS124  (CS554 ) 

GPSCS3  — >  GPS  -  CS017  (CS349) 

CS1CS2  — >  CS486  (CS529)  -  CS124  (CS554) 

CS1CS3  —  >  CS486  (CS529)  -  CS017  (CS349) 

THRESH-0.5D0 
INFILE  -  '  * 

OUTFILE  -  '  ' 

WRITE(*,*) 'PLEASE  ENTER  INPUT  USNO  (RDAS)  RAW  DATA  FILENAME 
READ( * , 101) INFILE 
FORMAT (A15) 

OPEN  ( 15 ,  FILE-INFILE ,  ACCESS- '  SEQUENTIAL '  , 

FORM- '  FORMATTED ' ,  STATUS-  •  OLD '  , ERR-999 ) 

WRITE (*,*) 'PLEASE  ENTER  OUTPUT  DATA  FILENAME 
READ ( * , 102 ) OUTFILE 
FORMAT (A15) 

OPEN ( 16 , FILE-OUTFILE , ACCESS- ' SEQUENTIAL '  , 

FORM-  '  FORMATTED ' ,  STATUS- '  UNKNOWN  ' ) 

READ (15, 9000, END-100 ) INDAT 

IF  (INDAT(l) .EQ. . AND . INDAT ( 2 ) .EQ. )  THEN 
READ (15,*) MJDAY 
READ (15,*) SYSTIME 
IF  (MOD (SYSTIME , 100) . NE. 0 )  GOTO  10 
MJDAY-DINT (MJDAY) 

WRITE (16,105) MJDAY , SYSTIME 
FORMAT ( IX , F10 . 1 , 3X, 17) 

READ ( 15 , * ) GPSCS2 
READ (15, *)GPSCS3 
IF  (GPSCS2 .GE. THRESH)  THEN 

GPSCS2— 10 .0D0**6*(GPSCS2-1.0D0) 

ELSE 

GPSCS2— 10 . 0D0**6*GPSCS2 
END  IF 

IF  ( GPSCS 3 .GE. THRESH)  THEN 

GPSCS3-10 . 0D0**6* (GPSCS3-1 . 0D0) 

ELSE 

GPSCS 3  =  10 . 0D0*  *  S*GPSCS3 
END  IF 

WRITE ( 16, *)GPSCS2,GPSCS3 

READ (1 5, *) DUMMY 

READ ( 1 5, *) DUMMY 

READ ( 15 , * ) CS1CS2 

READ ( 15 , * ) CS1CS3 

IF  (CS1CS2.GE. THRESH)  THEN 

CS ICS 2— 10 . 0D0**6* (CS1CS2-1 . 0D0 ) 

ELSE 


FILEPROC. FOR 


APRIL  21,  1990 


PAGE:  2 


CS1CS2-10. 0D0**6*CS1CS2 
END  IF 

IF  (CS ICS 3. GE. THRESH)  THEN 

CS1CS3“10. 0D0**6* (CS1CS3-1 . 0D0 ) 

ELSE 

CS1CS3»10 . 0D0**6*CS1CS3 
END  IF 

WRITE ( 16 ,*) CS 1CS2 ,  CS  ICS 3 
END  IF 
GOTO  10 

9000  FORMAT ( 8 0A1) 

100  WRITE (*, 9010)  OUTFILE 

9010  FORMAT ( *  STOPPED  PROCESSING  ',A15,'  FILE') 

CLOSE (15) 

CLOSE (16) 

STOP 

999  WRITE ( * ,  9020)  INFILE 

9020  FORMAT ( '  ERROR  OPENING  USNO  DATA  INPUT  FILE  ('^15,')') 
WRITE (*,*)  'PLEASE  CHECK  FILENAME  AND  TRY  AGAIN' 

END 


nono  n  o  o  o  oo  o  ooooooo 


GPS SIM. FOR 
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CCCSNOFLOATCALLS 

PROGRAM  GPSSIM 
WRITTEN  LBS  3/5/89 


THIS  PROGRAM  IS  A  COVARIANCE  SIMULATION  OF  CESIUM  CLOCK 
STABILITY  STUDY 


REAL* 8  TIME, DT,  P ( 10 , 10) , H ( 10) , R, F ( 10 , 10) , PHI ( 10 , 10) , 
Q(10 , 10) , QC (10,10) ,K(10) ,X(10) ,W1(10,10)  , 

W2 ( 10 , 10)  ,WT(10, 10) , PHASE (3) , FREQ (3) ,CALTIME(3) , 
MEASTIME ( 3 ) ,PER10D(3) ,FJUMP(3) , SIGGPS , SIGCLK( 3 ) 

CHARACTER* 5  CC 
CHARACTER* 15  FILE1 , FILE2 


COMMON/ CLKMOT/ PHASE , FREQ 
COMMON/ CLKPRM/FJUMP, PERIOD 
COMMON/MEASIG/ SIGGPS , SIGCLK 


N-6 
LSW-0 
FILE1- '  ' 

FILE2* '  ' 

OPEN (6, FILE- 'PROC. OUT' .STATUS* 'UNKNOWN ’ ) 

WRITE (6,2) 

2  FORMAT ( IX, '***  START  OF  PROGRAM  GPSSIM***') 

WRITE (*,*) 'PLEASE  ENTER  RAW  MEASUREMENT  DATA  FILE' 
READ ( * , 11 ) FILE2 

WRITE (*,*)' PLEASE  ENTER  NAME  FOR  LOTUS  OUTPUT  FILE' 
READ(* , 11) FILE1 
11  FORMAT (A15) 

OPEN ( 10 , FILE* ' PARM  .IN', ACCESS* ' SEQUENTIAL ' , 

1  FORM- ' FORMATTED ' , STATUS* ’ OLD ' , ERR-999) 

OPEN ( 20 , FILE-FILE1 , ACCESS- ' SEQUENTIAL '  , 

1  FORM- ' FORMATTED ', STATUS* 'UNKNOWN ' ) 

READ (10,*)  CC, (PHASE (I) , FREQ (I) ,1*1,3) 

WRITE (6  *)  'INITIAL  CLOCK  ERROR  ESTIMATES' 

WRITE (6, *)  (PHASE (I) , FREQ (I) ,1=1,3) 

DEFAULT  VALUES 

TIME-0. 

DT-0. 

DO  20  1*1,10 
H(I) =0. 0 
K( I) =0 . 0 
X(I)-0.0 
DO  10  J-l, 10 
P(I,J)«0. 

F ( I , J) — 0 . 

PHI ( T , J) -0 . 
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Q(I,J)-0. 

QC(I,J)-0. 

W1(I,J)«0. 

10  CONTINUE 
20  CONTINUE 
C 

c 

C  SET  INITIAL  CLOCK  ESTIMATES 

C 

DO  30  I-l,N/2 

X(2*I-1)  -  PHASE ( I ) 

X(2*I)  -  FREQ (I) 

30  CONTINUE 

C 

c 

READ ( 10 , *)  CC, (MEASTIME(I) ,1-1, N/2) 

WRITE ( 6 , * )  ' SAMPLING  FREQ  * 

WRITE ( 6 ,  * )  (MEASTIME(I) ,1-1, N/2) 

C 

READ(10, *)  CC, (P(I,I) ,1-1, N) 

C 

WRITE (6,*)  'INITIAL  COVARIANCE' 

C 

C  DO  40  1*1 , N/2 

c  P(2*I-1, 2*1-1)  -  P(2*I-1, 2*1-1)+ 

C  +  P(2*I , 2*1) 

c  40  CONTINUE 

WRITE (6, *) (P(I,I) ,1-1, N) 

C 

c 

C  READ  CLOCK  ERROR  MODEL  PARAMETERS: 

C  FREQUENCY  JUMP  MAGNITUDE 

C  JUMP  INTERVAL 

C 

READ (10,*)  CC, (FJUMP(I)  ,PERIOD(I) ,1-1,  N/2) 

CC  IF  NECESSARY  FREQUENCY  MAGNITUDE  JUMP  IS  ADJUSTED 
CC  TO  FIT  THE  VALUE  OF  0.0055  uSEC/DAY  (RMS),  BY  MAKING  A 
CC  FJUMP  *  0.03  USEC/DAY  IN  PARM.IN  FILE 

C 

C  CALCULATE  PROCESS  NOISE 

C 

C 

READ ( 10 ,  * )  CC,  (Q  (I ,  I)  ,1-1,  N’, 

DO  50  1-1 , N/2 

Q (2*1 , 2 *1 )  -  FJUMP(I) **2/(PERI0D(I)/MEASTIME(I) ) 
50  CONTINUE 
C 

WRITE (6,*)  'PROCESS  NOISE' 

C 

WRITE (6, *)  (Q ( I , I ) .1-1, N) 

C 

C 

READ ( 10 , * )  CC , SIGGPS , (SIGCLK(I) ,1-1, N/2) 

WRITE (6,*)  'MEASUREMENT  VARIANCE' 

C 

C  CALCULATE  MEASUREMENT  VARIANCE 

C 


SIGGPS  -  SIGGPS**2 


0 
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DO  60  I-l,N/2 

SIGCLK(I)  »  SIGCLK(I) **2 
60  CONTINUE 
C 

WRITE ( 6 , * )  SIGGPS, (SIGCLK(I) ,1-1, N/2) 

C 

C  RUN  KALMAN  FILTER 

C 

C 

CALL  KALMAN (TIME, DT, P,H,R, K, F,  PHI, QC,Q,N,W1,W2 , WT,N,X, FILE2) 
CLOSE (6) 

CLOSE (10) 

CLOSE (20) 

C 

C*****************END  MAIN  GPSSIM*********************** 

c 

STOP 

998  WRITE (6,*)  'ERROR  OPENING  RAW  TIMING  DATA  FILE' 

999  WRITE (6, *)  'ERROR  OPENING  PARM.IN  FILE' 

END 

C 

C 

SUBROUTINE  GETHR (MJDAY , TIME ,H,R,N, MOD2  ,  MOD3  ,  PARMS  ,  MEAZ  ,  W4  , 

1  FILE2) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  H,R  MATRICES  FOR  AN  UPDATE 
C  PLUS  IT  READS  THE  DATA  FILE  FOR  TIMING  MEASUREMENTS 

C 

CHARACTER* 15  FILE2 

REAL* 8  H ( 10) , R , MEAZ , TIME , MJDAY , Z ( 5 ) ,W4 (10) 

C 

C 

INTEGER  MOD2 , MOD3 , I 1 , I 2 , K1 , K2 , PARMS (80) 

INTEGER  SYSTIME 

C 

WRITE (6,*)  'NOW  IN  GETHR' 

WRITE ( 6 , ’ (A4) ' )  MOD2 
C 

C  INITIALIZE  H  TO  ZERO 
C 

DO  10  1-1,10 
H ( I ) *0 . 0 
10  CONTINUE 
C 
C 

IF (MOD2 . EQ . ' GPS1 ' )  CALL  HGPS1(H,R) 

IF (MOD2 .  EQ .  ' GPS2 '  )  CALL  HGPS2(H,R) 

IF (MOD2 . EQ. ' GPS3 ' )  CALL  HGPS3(H,R) 

IF (MOD2 . EQ . ' CLK1 ' )  CALL  HCLK1 (H , R) 

IF (MOD2 . EQ . ' CLK2 ' )  CALL  HCLK2 (H , R) 

C 

CALL  KFGETZ (MJDAY , TIME , SYSTIME , Z , W4 , FILE2 ) 

C 

IF (MOD2 . EQ . ' GPS1 ' )  MEAZ=Z(1) 

IF (MOD2 . EQ. ' GPS2 ' )  MEAZ=Z(2) 

IF(MOD2.EQ. 'GPS3 ' )  MEAZ=Z(3) 

IF (MOD2 . EQ . ' CLK1 ' )  MEAZ-Z(4) 

IF  (MOD2 . EQ . ' CLK2 ' )  MEAZ-Z(5) 
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C 

WRITE (6,*)  '  MEAS  DAY :  ' , MJDAY , 1  MEAS  TIME- SYSTIME , 
1  '  MEAS*  ' ,  MEAZ 

RETURN 
C 

END 

*******************END  SUBROUTINE  GETHR****************** 
SUBROUTINE  KFGETZ (MJDAY , TIME , SYSTIME , Z , W4 , FILE2 ) 
SUBROUTINE  TO  READ  AND  INTERPRET  THE  EVENT  TABLE. 


INTEGER  SYSTIME , DUM1 , LCOUNT 
R£AL*4  ZTIME , DUM2 , DUM3 

REAL* 8  Z(5) , TIME, SAVDT, THRESH, MJDAY, W4 (10) 
REAL* 8  GPSCS2 , GPSCS3 , CS1CS2 , CS ICS 3 
CHARACTER* 15  FILE2 

THRESH  »  0.5D0 

WRITE ( 6 #  * )  1  IN  KFGETZ  SUBROUTINE' 

OPEN ( 22 , FILE*FILE2 , ACCESS* ' SEQUENTIAL ' , 

1  FORM- ' FORMATTED ' , STATUS* ' OLD ' , ERR-999 ) 

CC  LCOUNT  *  IDINT(TIME) *24 

LCOUNT  -  TIME*24 . ODO 
DO  10  1*1, LCOUNT 

READ ( 2 2 , * ) DUM2 , DUM1 
READ ( 2 2 , * ) DUM2 , DUM2 
READ ( 2 2 , * ) DUM3 , DUM3 
10  CONTINUE 

READ ( 22 , * ) MJDAY , SYSTIME 
READ (22,*) GPSCS 2 , GPSCS  3 
READ ( 22 , *) CS1CS2 , CS1CS3 

2(1)  *  GPSCS2  -  CS1CS2 
Z  (2)  *  GPSCS 2 
2(3)  *  GPSCS 3 
2(4)  *  CS1CS2 
Z ( 5)  «  CS1CS3 

DO  20  1-1,5 
W4 (I) —Z (I) 

20  CONTINUE 


CLOSE (22) 

RETURN 

******************** END  SUBROUTINE  KFGETZ************************** 

999  WRITE (99 , * ) '  WARNING  ::  ERROR  OPENING  NORWAY . DAT  FILE' 

RETURN 

END 

C 

SUBROUTINE  FRTO (PARMS , II , 12 , K1 , K2 ) 
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C 

C  THIS  SUBROUTINE  FINDS  THE  'FROM'  AND  'TO'  VARIABLES  IN  UPDT 

C 

C 

CHARACTER  CPARMS (80) , FROM, TO 
INTEGER  II , 12 , K1 , K2 , ISTART , ILOC , PARMS (80) 

C 

FROM- ' F ' 

TO- ' T ' 

ISTART— 1 
ILOC— 0 
C 

CFIND  FIRST  FROM  INTEGER  AND  PUT  IN  II  VARIABLE 

C 

C 

C  THIS  CONVERTS  PARMS  AS  INTEGER  TO  CHARACTER  (I  HOPE) 

C 

DO  5  1-1,80 

CPARMS ( I )  -  CHAR ( PARMS ( I ) ) 

5  CONTINUE 
C 

c 

c 

CALL  SEARCH (FROM, CPARMS, ISTART, II, ILOC) 

IF(ILOC.EQ.O)  GO  TO  10 
C 

CNOW  GET  12 
C 

FROM- ' , ' 

I START- 1 
ILOC-O 

CALL  SEARCH (FROM, CPARMS, ISTART, 12, ILOC) 

C 

C  GET  TO  VARIABLE 
C 

I START- 1 
ILOC-O 

CALL  SEARCH (TO, CPARMS, ISTART, Kl, ILOC) 

C 

C  GET  NEXT  ANTENNA 
C 

TO* ' + ' 

ISTART— 1 
ILOC-O 

CALL  SEARCH (TO, CPARMS, ISTART, K2, ILOC) 

C 

WRITE (6 , * )  'SAT.  AND  ANT.  FOR  THIS  MEAS . ' 

WRITE (6,'(4(1X,I2))')  11,12, Kl , K2 
C 

RETURN 

C 

10  WRITE (6, 100) 

100  FORMAT (IX,'  PARM  CARD  FORMAT  ERROR  IN  UPDT') 

RETURN 

C 

END 

C 

c****************end  of  frto***************************** 
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C 

SUBROUTINE  SEARCH ( CHAR , STRING , ISTART , I VAR , ILOC) 


THIS  ROUTINE  SEARCHES  'STRING'  FOR  'CHAR'  STARTING  AT 
ISTART.  IT  THEN  SEARCHES  FOR  THE  NEXT  EQUAL  SIGN  ANS 
SKIPS  UP  TO  FIVE  CHARACTERS  TO  FIND  THE  INTEGER  AND  PUT  IT 
IN  I  VAR.  IF  CHAR  OF  EQ  NOT  FOUND,  THEN  ILOOO 

CHARACTER* 1  STRING ( 8 0 ), EQ , CHAR 
INTEGER  ISTART, IVAR, ILOC, INDEX 


EQ* • » • 

FIND  ’CHAR’  IN  'STRING* 

CALL  FINDC( CHAR, STRING, ISTART, ILOC) 

IF ( ILOC. EQ . 0)  RETURN 

ISTART-ILOC 

IL0C*0 

FIND  EQ  SIGN 

CALL  FINDC(EQ, STRING, ISTART, ILOC) 

IF (ILOC. EQ. 0)  RETURN 

SEARCH  UP  TO  FIVE  SPACES  FOR  NONBLANK  CHAR 
DO  10  INDEX  *1,5 

IF (STRING* ILOC+INDEX) .NE. '  »)  GO  TO  20 
10  CONTINUE 

ILOC=0 

RETURN 

20  READ (STRING (ILOC+INDEX) , 100)  IVAR 
100  FORMAT (II) 

RETURN 

END 


SUBROUTINE  KALMAN (TIME , DT , P , H , R , K , F , PHI , QC , Q , N , W1 , W2 , WT , M , X , FI LE: 


THIS  SUBROUTINE  DOES  COVARIANCE  ANALYSIS  VIA  KALMAN  FILTERING 
IMPLICIT  R£AL*8  (A-H,P-Z) 


REAL*8  P(10,10) ,F(1Q,10) , PHI (10,10) ,H(10) ,QC(10,10) ,Q(10,10) 
REAL* 8  Wl(10,10) ,W2(10,10) ,WT(10,10) , X ( 10 ) ,W3(10) ,W4(10) 
REAL* 8  K  ( 10)  , NEWDT ,  EVTAB  (50,8)  , MJDAY  ,  nrF , TIME 
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INTEGER  EVENT,  P ARMS  (80)  , MODI , M0D2  , M0D3  ,  PRMTAB (50,80)  ,N 
REALM  SEVTAB(50,8)  , STIME , SSDT,  SNEWDT 
CHARACTER* 15  FILE 2 


SET  UP  EVENT  TABLE 

CALL  KFTABL(SEVTAB, PRMTAB) 

DO  3  1-50 
DO  3  J— 8 

EVTAB(I, J)  -  DBLE ( SEVTAB ( I , J ) ) 
3  CONTINUE 


SET  INITIAL  TIME  AND  DT 
TIME-0. 0D0 
DT-EVTAB (21,7) 

READ  CONTROL  CARD 
STIME  -  0.0 
SSDT  -  SEVTAB (21,7) 

1000  CALL  KFREAD( STIME, SSDT, SNEWDT, SEVTAB,  PRMTAB, EVENT, 

1  MODI , MOD2 , M003 , PARMS ) 

DT  -  DBLE (SSDT) 

NEWDT  -  DBLE (SNEWDT) 

END  OF  FILTER  RUN 

IF  ( EVENT. EQ. 'STOP' )  GOTO  9999 

PROPAGATE 

IF  (EVENT. EQ. 'PROP' )  CALL  KFPROP (MJDAY , TIKE , DT, NEWDT, MODI , 

1  MOD2  , MOD3  ,  N ,  P , H ,  F ,  PHI ,  QC ,  Q ,  PARMS  , W1 , W2 ,W3 , W4  ,  X,  FILE2) 

IF ( EVENT. EQ. 'PROP' )  STIME  *  SNGL(TIME) 

OR  UPDATE 

IF  ( EVENT. EQ.  'UPDT' ) CALL  KFUPDT ( MJDAY, TIME  MOul ,MOD2 , MOD3 , N , P ,  K , 

1  R, K, PARMS , W1 , W2 , W3 , W4 , X, FILE2 ) 

OR  PRINT 

IF  (EVENT. EQ. ' PRNT ' )  CALL  KFPRNT (TIME , MODI , MOD2 , MOD3 , P , F , PHI , H , R, 
1  K , QC , Q , N , PARMS , W1 , X ) 

OR  USER  DEFINED 

IF  (EVENT  .EQ.  'USER')  CALL  KFUSER (MJDAY , TIME , DT , MODI , MOD2 , MOD3  , 

1  P,  F,  PHI ,  Q,  QC,  H,  K,  R,  WT,  W1 ,  W2 ,  W3  ,  W4  ,  PARMS,  M,N,X) 

GOTO  1000 
9999  RETURN 
END 


******************£1^  OF  KALMAN  ROUTINE  **•+*+**+*******•***•**•** 
SUBROUTINE  ABADD ( A , B , M , N , C ) 
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SUBROUTINE  TO  DO  THE  MATRIX  ADDITION  C-A+B 
THE  MATRICES  A,B,C  NEED  NOT  BE  DISTINCT 

IMPLICIT  REAL* 8  (A-H.P-Z) 

DIMENSION  A(10 , 10) ,B(10,10) ,C(10,10) 

DO  10  J-1,N 
DO  10  I-1,M 

10  C(I,J)-A(I,J)+B(I,J) 

RETURN 

END 

SUBROUTINE  ABMULT (A, B, L,M,N , C , Wl) 


SUBROUTINE  TO  COMPUTE  THE  MATRIX  PRODUCT  C=A'B 

A, B , C  NEED  NOT  BE  DISTINCT.  Wl  CAN  NOT  BE  THE  SAME 

AS  A  OR  B 

IMPLICIT  R£AL*8  (A-H,P-Z) 

REAL*8  A ( 10 , 10) , B ( 10 , 10) ,C(10,10) ,W1(10,10) 

DO  20  J«1,N 
DO  20  I«1,L 
SUM-0.0 DO 
DO  10  K-l ,M 

10  SUM-SUM+A ( I , K) *B(K,  J) 

20  Wl (I , J) —SUM 

DO  30  J-1,N 
DO  30  1-1 , L 

30  C(I,J)«W1(I,J) 

RETURN 

END 

SUBROUTINE  KFERR (TIME , CARD) 

IMPLICIT  REAL* 8  (A-H,P-Z) 

SUBROUTINE  TO  STOP  THE  PROGRAM  IN  CASE  OF  A  SERIOUS  ERROR. 
INTEGER  CARD (8) 

WRITE(6, 100)  TIME, (CARD(I) ,1-1,8) 

100  FORMAT (///, '*****  IMPROPER  EVENT  AT  TIME-  ',G10.5, 

1  '  EVENT  CARD  IS:  ' ,/,5X,8A10) 

STOP 

END 

SUBROUTINE  KFHPHT (A, B , M, N, C, Wl ) 


SUBROUTINE  TO  COMPUTE  THE  MATRIX  PRODUCT  C=A*B*TRANS (A) 
WHERE  B  IS  SYMMETRIC  AND  THE  RESULT  IS  SYMMETRIC. 

THE  MATRICES  A,B,C  NEED  NOT  BE  DISTINCT. 

IMPLICIT  REAL* 8  (A-H,P-Z) 

R£AL*8  A (10, 10) ,B(10,10) ,C(10,10) ,W1(10,10) 
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200 


300 

100 

10 


DO  10  I-1,M 
DO  10  J-I,M 
TERM-0. 0D0 
DO  100  K— 1,N 
SUM-0 . 0D0 
DO  200  L-l , K 

SUM-SUM+A (I , L) *B(L,K) 
IF  (K.EQ.N)  GOTO  100 
KP1  -  K+l 
DO  300  L-KPl.N 

SUM-SUM+A ( I , L) *B (K, L) 
TERM-TERM+SUM* A ( 3 , K) 

HI (I , J) —TERM 


COPY  HI  INTO  RETURN  MATRIX  C 


DO  20  I— 1,M 
DO  20  J— I ,M 

C(I , J) —HI (I, J) 

20  C(J,I)«W1(I,J) 

RETURN 

END 

SUBROUTINE  KFPRNT ( TIME , MODI , MOD2 , MOD3 , P , F , PHI , H , R , K , QC , Q , N , PARMS , 
1  W1,X) 


SUBROUTINE  TO  PRINT  THE  DESIRED  MATRIX  IN  IT’S  DESIRED  FORM. 
IMPLICIT  REAL* 8  (A-H,P-Z) 


REAL* 8  P(10,10) ,F(10,10) , PHI (10, 10) ,QC(10, 10) ,Q(10, 10) 
REAL* 8  H ( 10) , W1 (10,10) 

REAL* 8  K(10) ,X(10) 

INTEGER  PARMS ( 80 ), MODI, MOD2 ,M0D3 


WRITE (6, 100)  TIME , MODI , MOD2 , (PARMS (J) ,J=1,80) 

100  FORMAT (//, '***  TIME-  ' ,G10. 3 , 2 (5X, A4 ) , IX, 80A1 ) 

DETERMINE  PRINT  TYPE 
IF  (MODI. EQ. 'FULL' )  GOTO  1000 
IF  (MODI . EQ . ' CORR ' )  GOTO  2000 
IF  (MODI . EQ . ' LTRI ' )  GOTO  3000 
IF  ( MODI . EQ . ' PRMS ' )  GOTO  4000 
IF  ( MODI . EQ . ' UDER ' )  GOTO  5000 
IF  (MOD1.EQ. 'DEBG' )  GOTO  6000 

UNDEFINED  PRINT  TYPE.  PRINT  WARNING  AND  RETURN. 

WRITE (6,200) 

200  FORMAT (//, ******  UNRECOGNIZED  PRINT  OPTION  SPECIFIED  **♦**’) 
GOTO  9999 

PRINT  THE  FULL  MATRIX 
1000  IF  (MOD2.EQ. ’ P ' )  CALL  PRINT(P,N,N) 

IF  (M0D2.EQ. 'F')  CALL  PRINT (F,N,N) 

IF  (M0D2.EQ. 'X' )  CALL  VPRINT(X,N) 
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IF  (MOD2.EQ. 'PHI')  CALL  PRINT (PHI, N,N) 

IF  (MOD2.EQ. 'Q')  CALL  PRINT (Q,N,N) 

IF  (MOD2.EQ. 'QC')  CALL  PRINT (QC,N,N) 

IF  (MOD2.EQ. * H * )  CALL  VPRINT(H,N) 

IF  (MOD2.EQ. * K* )  CALL  VPRINT(K,N) 

IF  (MOD2.EQ. 'R«)  WRITE(6 , 1100)  R 
1100  FORMAT (/,5X,G10. 3) 

GOTO  9999 
C 

C  PRINT  THE  CORRELATION  COEFFICIENTS 

2000  IF (MOD2 . EQ . ' P 1 )  GOTO  2002 

2001  WRITE(6,2101) 

2101  FORMAT (//, ******  IMPROPER  MATRIX  SPECIFIED  *****') 

GOTO  9999 

2002  DO  2022  I»1,N 
DO  2012  J»1,I 

2012  W1(I,J)«P(I,J)/DSQRT(P(I,I) *P(J,J) ) 

2022  WRITE (6,2122)  I , (W1 ( I , J) , J-l , I ) 

2122  FORMAT (// ,I5,3X,10(G10.3,1X)/(8X,10(G10.3,1X) ) ) 

GOTO  9999 
C 

C  PRINT  THE  LOWER  TRIANGULAR  PORTION  OF  THE  MATRIX 
3000  GOTO  9999 

C 

C  PRINT  THE  RMS  VALUES  OF  THE  P  MATRIX 
4000  IF  (MOD2.NE.'P  ’)  GOTO  2001 

DO  4010  I-1,N 

4010  W1(1,I)-SQRT(P(I,I)  ) 

WRITE  (*,4100)  (W1(1,I) ,I«1,N) 

4100  FORMAT (// ( 3X, 10 (IX , G10 .  3 )  )  ) 

GOTO  9999 

C 

C  PRINT  USER  DEFINED  MATRICES 

5000  GOTO  9999 

CC  CALL  UPRINT(TIME,MOD2 , MOD3 , PARMS , P, F, PHI , Q , QC, H , K, N , W1 ) 

C 

C  PRINT  DEBUG  MATRICES 

6000  GOTO  9999 
C 

9999  RETURN 

END 

SUBROUTINE  KFPROP  (MJDAY , TIME , DT , NEWDT , MODI , M0D2 , MOD3 , N , P , H , 
1  F,  PHI, QC,Q, PARMS , Wl, W2 ,W3,W4 ,X,FILE2) 

C 

C 

C  SUBROUTINE  TO  DO  KALMAN  FILTERING  COVARIANCE  PROPAGATION 

IMPLICIT  REAL* 8  (A-H,P-Z) 

C 

C 

REAL* 8  P(10, 10) ,F(10, 10)  ,  PHI (10, 10) , QC ( 10 , 10) , Q ( 10 , 10) 

REAL* 8  Wl(10,10) ,W2( 10,10) ,W3( 10) ,H(10) ,MEAZ(5) ,W4(10) 

REAL* 8  X ( 10 ) , MJDAY , TIME , Z , HX 

INTEGER  MODI, MOD2,MOD3, PARMS (80) ,SYSTIME,N 

REAL* 8  NEWDT , OLDDT , DT 

CHARACTER* 15  FILE2 

DATA  OLDDT/ 0.0/ 

C 

c 
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WRITE (6 ,  *)  '  HOW  IN - KFPROP - 1 

C 

C  CHECK  FOR  TIME  VARYING  SYSTEM 
C 

CCCC  IF  (MODI. EQ. 'VARY* )  CALL  GETFQC  (TIME,  DT,  F,  QC, N ,MOD2  ,MOD3  ,  PARMS) 

C 

C  COMPUTE  NEW  PHI  AND  Q  IF  NEW  DT  OR  NEW  F  OR  NEW  QC 
C 

IF  ( (OLDDT  .NE.  DT)  .OR.  (MODI  .EQ.  'VARY')) 

1CALL  PHICOM(DT, F , QC, Q, N, PHI , W1 , W2 ) 

C 

C  PROPAGATE  THE  STATE, IF  REQUESTED 

C 

IF  (MOD2.NE. 'STAT')  GOTO  30 
DO  10  1*1, N 

W1(I,1)«0.0 
DO  20  J-1,N 

W1 ( I , 1) *W1 (1,1) +PHI (I , J) *X(J) 

20  CONTINUE 

10  CONTINUE 

DO  25  1*1  ,N 
X(I)-W1(I,1) 

25  CONTINUE 
C 

C  INSERTED  BY  ATO  4/15/85  TO  SIMULATE  DRIVING  TERM 
cccc  IF (MOD3 . EQ . ' DRVN ' )  CALL  GETU(X,N,DT) 

C 

C 

C  PROPAGATE  P  MATRIX 

C 

30  CALL  KFHPHT ( PHI , P , N , N , P , W1 ) 

CALL  ABADD(P,Q,N,N, P) 

C 

C  PROPAGATE  CLOCK  READOUT 
C 

CCCCC  CALL  RHOCOM ( TIME, DT, MOD2 ,MOD3 , PARMS) 

C 

C  CALCULATE  MEASUREMENT  RESIDUALS  AT  PROPAGATE  TIME 

C 

C 

c 

CALL  KFGETZ (MJDAY , TIME , SYSTIME , MEAZ , W4 , FILE2 ) 

C 

DO  47  1*1,5 
DO  45  J*1,N 
H ( J)  *  0.0D0 
4  5  CONTINUE 

IF  (I.EQ.l)  THEN 
H ( 1) *-l . 0D0 
ELSE 

IF  (I.EQ.2)  THEN 
H(3) *-l. 0D0 
ELSE 

IF  (I.EQ.3)  THEN 
H(5)— l.ODO 
ELSE 

IF  (I.EQ.4)  THEN 
H ( 1 ) *1 . 0D0 
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H(3)»- 1.0D0 

ELSE 

IF  (I.EQ.5)  THEN 
H(1)«1.0D0 
H(5)— 1.0D0 
END  IF 
END  IF 
END  IF 
END  IF 
END  IF 
HX-0.0D0 
DO  46  J-1,N 

HX-HX+H(J)*X(J) 

4  6  CONTINUE 
C 

W3  (I)  -  MEAZ(I)  -  HX 
47  CONTINUE 

WRITE ( 6 ,  * )  '  MEAS  DAY : ' , MJDAY , '  MEAS  TIME- ' , SYSTIME , 

1  '  MEAS-  ' ,  MEAZ 

WRITE (6,*)'  STATE  VECTOR  X(I) :  ',X 
WRITE (6, *) '  RESIDUALS:  ' ,W3 

UPDATE  TIME  AND  DT 

TIME  -  TIME  +  DT 
OLDDT  -  DT 
DT  -  NEWDT 
RETURN 
END 

SUBROUTINE  KFUPDT (MJDAY , TIME  r MODI , MOD2 , MOD3 , N , P , H , R , K , PARMS , 
1  W1 , W2 , W3 , W4 , X , FILE2 ) 


SUBROUTINE  TO  PERFORM  THE  KALMAN  FILTER  COVARIANCE 
MATRIX  UPDATE. 

W3  vector  contains  the  residuals  for  measurements  used 
in  this  run  of  a  filter 

IMPLICIT  REAL* 8  (A-H,P-Z) 


REAL*8  P(10,10) ,H(10) ,X(10) ,W1(10,10) ,W2(10,10) 
INTEGER  MODI , MOD2 , MOD3 , PARMS ( 8  0 )  ,  M ,  N 
REAL* 8  K(10) , MJDAY, TIME, W3 (10) ,W4(10) ,Z,HX 
CHARACTER* 15  FILE2 


WRITE (6,*)'  NOW  IN  -  KFUPDT  - ' 

GET  NEW  H  AND  R  IF  NEEDED. 

IF  (MODI. EQ. ’VARY' ) 

1  CALL  GETHR  (MJDAY,  TIME,  H,R,N,MOD2,MOD3,  PARMS  ,Z,W4,FILE2) 
COMPUTE  THE  KALMAN  GAINS  MATRIX,  K. 

HPHT-0.0D0 
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CO  20  J-1,N 
SUM-0. 0D0 
DO  10  1-1, N 

10  SUM— SUM+H ( I ) *P ( I >  J ) 

20  HPHT-HPHT+SUM*H(J) 

C 

HPHT— HPHT+R 
DO  40  1*1, N 
SUM-0. 0D0 
DO  30  J-1,N 

30  SUM— SUM+P(I , J) *H(J) 

40  K ( I ) -SUM/HPHT 

IF  (MOD3.NE. 'STAT' )  GOTO  49 

C 

C  UPDATE  STATE  VECTOR 

C 

HX-0.0D0 
DO  47  I-l , H 

HX-HX+H ( I ) #X ( I ) 

47  CONTINUE 
C 

DO  48  I— 1 , N 

X(I) — X (I) +K( I) * ( Z-HX) 

48  CONTINUE 
CC 

C  CALCULATE  RESIDUAL 

C 

HX-0.0D0 
DO  46  I— 1,N 

HX-HX+H (I) *X(I) 

46  CONTINUE 

C 


IF 

(MOD2.EQ. 'GPS1 ' ) 

W3(l) 

at 

Z  -  HX 

IF 

(MOD2.EQ. ' GPS2 ' ) 

W3(2) 

* 

Z  -  HX 

IF 

(MOD2.EQ. ' GPS 3 ' ) 

W3(3) 

ss 

Z  -  HX 

IF 

(MOD2.EQ. ' CLK1 ’ ) 

W3(4) 

* 

Z  -  HX 

IF 

(MOD2.EQ. ' CLK2 ' ) 

W3(5) 

Z  -  HX 

49  CONTINUE 
C 

C  UPDATE  THE  COVARIANCE  MATRIX,  P. 

C 

DO  60  I-l , N 
DO  50  J— 1 , N 

50  W1(I,J)—  K(I)*H(J) 

60  HI ( I , I) *W1 ( I , I ) +1 . ODO 

C 

CALL  KFHPHT (W1 , P, N , N , P, W2 ) 

DO  70  I-l , N 
DO  70  J— 1 , N 

70  W1(I,J)-K(I)*K(J)*R 

CALL  ABADD(P,W1,N,N,P) 

RETURN 

END 

SUBROUTINE  PHICOM ( DT , F , QC , Q , N , PHI , W1 , W2 ) 
C 
C 
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C  SUBROUTINE  TO  COMPUTE  PHI  *  EXP  (F*DT) 

IMPLICIT  R£AL*8(A-H,P-Z) 

C 

C 

REAL* 8  F(10,10) , QC ( 10 , 10) , Q ( 10 , 10 ) , PHI (10, 10) 
REAL* 8  Wl(10,10) ,W2(10,10) 

REAL* 8  L1NORM , K , KFACT 
DATA  DELTA , KCOUNT/ 1 . 0D- 5 , 2  0/ 

C 

C  COMPUTE  THE  LI  NORM  OF  F*DT 
C 

L1NORM-O.ODO 

M-0 

C 

C  DO  20  J-1,N 

C  TEMP  -  0.0D0 

C  DO  10  I-1,N 

C  10  TEMP-TEMP+DABS(F(I, J) ) 

C  20  IF  (TEMP . GT . L1NORM)  L1NORM-TEMP 

C  L1N0RM-L1N0RM*  DT 

CC 

C  CHECK  TO  SEE  IF  NORM  IS  ZERO 
C  IF  (L1NORM  .LT.  1.0D-10)  GOTO  25 

C 

C  SCALE  F*DT 

C 

C  RLOG  *  (DLOGIO (L1NORM) ) / (DLOGIO ( 2 . 0D0) ) 

C  IF  (L1NORM.LT. 1)  M-0 

C  IF  (L1NORM.GE.1  .AND.  L1NORM.LT.2)  M«1 

C  IF  (L1NORM.GE.2)  M-I DINT (RLOG) +1 

C  25  TEMP— DT/ ( 2 . 0D0**M) 

C  DO  40  J— 1 ,N 

C  DO  30  I— 1 , N 

C  PHI (I , J) — F(I , J) *TEMP 

C  30  W1(I,J)-PHI(I,J) 

C  40  PHI ( J , J) -PHI ( J , J) +1 . 0D0 

C  L1NORM-L1NORM/ ( 2 . 0D0**M) 

CC 

C  COMPUTE  EXP(Wl)  BY  TAYLOR 

C 

C  K-1.0D0 

C  KFACT- 1 . 0D0 

C  1000  K-K+1.0D0 

C  KFACT-KFACT*K 

C  CALL  ABMULT ( F , W1 , N , N , N , W1 , W2 ) 

C  DO  60  J-l , N 

C  DO  60  I— 1 , N 

C  W1(I,J)«W1(I, J) *TEMP 

C  60  PHI ( I , J) —PHI (I,J)+W1(I,J) /KFACT 

C  REM- ( L1NORM** (K+l . 0D0 ) ) / 

C  1  (KFACT* (K+l . 0D0) * ( 1 . ODO-L1NORM/ ( K+2 . 0D0 ) ) ) 

C  IF  (K.GT. KCOUNT)  GOTO  9000 

C  IF  (DABS (REM) .GT. DELTA)  GOTO  1000 

CC 

C  CONVERGENCE  HAS  OCCURRED 
C 

c 
c 
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DO  60  J-1,N 

PHI(J,J)-1.0D0 
60  CONTINUE 

DO  70  J-1,N,2 
PHI(J,J+1)-DT 
70  CONTINUE 

C 

C  COMPUTE  Q  OVER  SUBINTERVAL  DT/(2**M) 

C 

TEMP-DT/ (2 . 0D0**M) 

C 

CCCC  CALL  QC0MP(TEMP,F,PHI,N,QC,Q,W1,W2,M) 

C 

C 

C  WRITE  WARNING  INFO  INTO  FILE  99. 

CCC  WRITE (99, 17 9)  M,K 

CCC  199  FORMAT (/ '  ***  PHICOM  **•  M-',I3,'  #  OF  TERMS  »',F5.0) 

C 

C 

C  IF  M  «  0  THEN  RETURN 

C  IF  M  >  0  POWER  UP  PHI  2**M  TIMES 

C 

IF  (M.EQ.O)  GOTO  9999 
C 
C 

C  POWER  UP  PHI  2**M  TIMES 

C 

DO  510  I-1,M 

510  CALL  ABMULT ( PHI , PHI , N , N , N , PHI , W2 ) 

C 

C  POWER  UP  COMPLETE 

9999  RETURN 

C 

C  CONVERGENCE  FAILURE  BLOCK.  STOP  PROGRAM. 

C 

9000  WRITE (6,9100) 

9100  FORMAT (///, 4 1H****  CONVERGENCE  FAILED  IN  PHICOMP  ****, 

1  // , 16H**  F  MATRIX  **) 

CALL  PRINT (F,N,N) 

STOP 

C 

END 

SUBROUTINE  PRINT (A, M, N) 

C 

C  SUBROUTINE  TO  PRINT  AN  M  BY  N  MATRIX,  A. 

C 

REAL*8  A ( 10 , 10) 

DO  10  I-1,M 

WRITE ( 6 , 100 )  I, (A(I, J) ,J»1,N) 

10  CONTINUE 

100  FORMAT (// ,I5,3X,10(G10.3,1X)  ,/ 

1  (8X,10(G10.3,1X) ) ) 

RETURN 

END 

SUBROUTINE  QCOMP ( DT , F , PHI , N , QC , Q , W1 , W2 , M ) 

C 

C 

C  SUBROUTINE  TO  COMPUTE  THE  DISCRETE  FORM  OF  THE 
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C  PROCESS  NOISE,  Q. 

C 

IMPLICIT  REAL*8  (A-H,P-Z) 

C 

REAL* 8  PHI (10, 10) ,F(10,10) ,QC(10,10) ,Q(10,10) 

REAL*8  Wl(10,10) ,W2(10,10) 

LOGICAL  FLAG 
C 

DATA  RELERR, ABSERR/1 . OD-6 , 1 . 0D-20/ 

C 

DO  50  J-1,N 
DO  50  I-1,N 

W1(I,J)-DT*QC(I,J) 

50  Q(I,J)-0.0 

C 

C  EXPAND  Q  OVER  INTERVAL  DT 

C 

DO  1000  ITER-2,25 
C 

C  ADD  EXPANSION  TERM  W1  TO  Q 

C  COMPUTE  W2-DT*F*W1 

C 

DO  100  I— 1 , N 
DO  100  J— 1 , N 

Q(I, J)-Q(I , J)+W1(I, J) 

SUM-0 . 0D0 
DO  90  K«1,N 

90  SUM— SUM+F ( I , K) *W1 ( K , J ) 

100  W2 (I , J) — SUM*DT 

C 

C  COMPUTE  NEXT  EXPANSION  TERM,  W1 

C 

SCALE- l.ODO/ ITER 
DO  200  J— 1 , N 

DO  200  I— 1 , N 

VALUE«(W2(I,J)+W2(J,I) ) *SCALE 
W1 ( I , J ) —VALUE 
200  W1 ( J, I) —VALUE 

C 

C  ZERO  OUT  NEGLIGIBLE  TERMS  IN  EXPANSION 
C 

FLAG-. FALSE. 

DO  300  J— 1 , N 
DO  300  I— 1 , N 

IF (DABS (W1 (I , J) ) .  LT. (ABSEER+RELER*DABS (Q ( I , J) ) ) )  GOTO  250 
FLAG- . TRUE . 

GOTO  300 

250  W1(I,J)«0.0D0 

3  00  CONTINUE 

C 

C  CONVERGENCE  CRITERION  IS  W1=0 
C 

IF(. NOT. FLAG)  GOTO  1100 
1000  CONTINUE 
1100  CONTINUE 
C 

C  IF  M  -  0  THEN  RETURN  ELSE  PROPAGATE  Q  ACROSS  INTERVAL  DT. 

C 
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IF  (M  .EQ.  0)  GOTO  1400 
MAX  -  2**M 
00  1200  J«1,N 
DO  1200  I-1,N 
1200  W1(I,J)-Q(I,J) 

DO  1300  1*1 , MAX 

CALL  KFHPHT ( PHI , Q , N , N , Q ,  W2 ) 

1300  CALL  ABADD(Q,W1,N,N,Q) 

1400  CONTINUE 

WRITE  WARNING  INFO  TO  FILE  99. 

CC  WRITE (99, 199)  ITER 

CC  199  FORMAT (/ '  ***  QCOMP  ***  ITER-', 15) 

RETURN 

END 

SUBROUTINE  VPRINT(V,N) 

SUBROUTINE  TO  PRINT  AN  N  VECTOR,  V. 

IMPLICIT  REAL* 8  (A-H,P-Z) 

REAL* 8  V(10) 

WRITE(6, 100)  (V(J) , J-1,N) 

100  FORMAT (// , ( 10 ( IX , G10 . 3 ) ) ) 

RETURN 

END 

FUNCTION  TRACE (A, M) 

FUNCTION  TO  COMPUTE  THE  TRACE  OF  THE  M  BY  M  MATRIX  A. 

IMPLICIT  REAL* 8  (A-H,P-Z) 

REAL*8  A(M,M) 

SUM  -  0.0 
DO  10  I— 1  ,M 
10  SUM  -  SUM  +  A(I, I) 

TRACE  -  SUM 

RETURN 

END 

SUBROUTINE  KFREAD ( TIME , DT , NEWDT , EVTAB , PRMTAB , EVENT , 

1  MODI , MOD2 , MOD3 , PARMS ) 


SUBROUTINE  TO  READ  AND  INTERPRET  THE  EVENT  TABLE. 


INTEGER  EVENT, MODI ,MOD2 ,MOD3 , PRMTAB (5 0,80) , PARMS (80) 
REAL  NEWDT, EVTAB (5 0,8) 

REAL  REVENT , RMOD1 ,  RMOD2 , RMOD3 

EQUIVALENCE  (I EVENT, REVENT) , ( IMOD1 , RMOD1) , ( IMOD2 , RMOD2 ) 
1  (IMOD3 , RMOD3 ) 


DT02  -  DT/2.0 
DO  100  INDEX-1,25 

IF  ( ABS ( TIME-EVTAB ( INDEX , 5 ) )  .GT.  DT02 )  GOTO  100 
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IF  (((TIME  .LT.  ( EVTAB  ( INDEX ,  6 )  +DT02 )  )  .OR. 

1  (EVTAB (INDEX, 6)  .EQ.  0.0))  .AND.  ( EVTAB ( INDEX , 1 )  .NE.  0.0) 

2  GOTO  200 

EVTAB ( INDEX , 5 ) -1 . 0E3 0 
100  CONTINUE 

NO  MATCHING  TIME,  TERMINATE  RUN 

EVENT  ■  'STOP1 
GOTO  9999 

TIME  MATCH,  SET  EVENT  AND  MODS 

200  REVENT  -  EVTAB ( INDEX , 1 ) 

RMOD1  -  EVTAB (INDEX, 2) 

RMOD2  -  EVTAB (INDEX, 3) 

RMOD3  -  EVTAB (INDEX,  4) 

COPY  EVENT  AND  MODIFIERS  INTO  PARAMETERS  FOR  RETURN 
EVENT  ■  I EVENT 
MODI  «  IMOD1 
MOD2  -  IMOD2 
MOD3  -  IMOD3 

SET  PARMS,  0  IN  EVTAB ( INDEX , 8 )  IMPLIES  NO  FARMS 

IPARM  -  EVTAB (INDEX, 8) 

IF  (IPARM  .EQ.  0)  GOTO  350 
DO  300  1-1,80 

PARMS ( I ) -PRMTAB ( I PARM , I ) 

300  CONTINUE 
GOTO  450 

350  DO  400  1-1,80 
PARMS ( I ) — •  ' 

400  CONTINUE 
450  CONTINUE 

UPDATE  EVENT  NEXT  TIME 

EVTAB (INDEX, 5)  -  EVTAB (INDEX, 5)  +  EVTAB ( INDEX , 7 ) 

IF  ( EVTAB ( INDEX , 7 ) . EQ . 0 . 0 )  EVTAB (INDEX, 5) -EVTAB (INDEX, 5) +  DT 
IF  (EVTAB (INDEX, 6)  .EQ.  0)  EVTAB ( INDEX, 5) -1 . 0E30 

IF  EVENT  IS  PROP  THEN  CHECK  FOR  NEW  DT 

IF  (EVENT  .EQ.  'PROP')  GOTO  500 

IF  EVENT  IS  AT  LAST  TIME  SET  NEXT  TIME  TO  INFINITY 

IF  (EVTAB (INDEX, 5)  .GT.  (EVTAB ( INDEX , 6 ) +DT02 ) ) 

1  EVTAB (INDEX, 5 )-1.0E30 
GOTO  9999 

PR PRO PAG ATE  EVENT  DT  CHECK 
500  CONTINUE 

IF  NOT  AT  THE  END  OF  A  PROPAGATE  SEQUENCE  THEN  RETURN 
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NEWDT  -  DT 

IF  (ABS(EVTAB( INDEX,  5)  -EVTAB(INDEX,  6)  )  .GT.  DT02)  GOTO  9999 


FIND  NEXT  PROPAGATE  SEQUENCE  IF  IT  EXISTS  AND  SET  NEW  DT 
DO  600  J-21,25 

IF  (ABS  (EVTAB(  INDEX,  5)  -EVTAB  ( J,  5)  )  .LT.  DT02  .AND. 
1  INDEX  .NE.  J)  GOTO  700 

600  CONTINUE 

NO  NEW  PROPAGATE  SEQUENCE 
NEWDT  -  DT 

EVTAB (INDEX, 5)  *  1.0E30 
GOTO  9999 

NEW  PROPAGATE  SEQUENCE.  SET  NEW  DT. 

700  NEWDT  -  EVTAB(J,7) 

EVTAB (INDEX, 5)  »  1.0E30 

999  CONTINUE 
RETURN 
END 


SUBROUTINE  HGPS1(H,R) 

THIS  SUBROUTINE  COMPUTES  THE  H  MATRIX  FOR  GPS  CLOCK  MEAS. 
GPS  -  CS486 

REAL* 8  H(10) , R, SIGGPS , SIGCLK (3 } 

COMMON/MEASIG/SIGGPS , SIGCLK 
COMPUTE  THE  NON-ZERO  ENTRIES  OF  H 

WRITE (6,*)  •  NOW  IN  HGPS1  SUBROUTINE' 

H ( 1) “-1 . 0D0 

MEAS.  NOISE 

R*DSQRT (SIGGPS **2+SIGCLK ( 1 ) **2) 

RETURN 


***************** END  SUBROUTINE  HGPS1************************* 


END 


SUBROUTINE  HGPS2(H,R) 

THIS  SUBROUTINE  COMPUTES  THE  H  MATRIX  FOR  GPS  CLOCK  MEAS. 
GPS  -  CS 12 4 
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REAL*8  H(10) , R, SIGGPS , SIGCLK ( 3 ) 

COMMON/MEASIG/SIGGPS , SIGCLK 
COMPUTE  THE  HOH-ZERO  ENTRIES  OF  H 

WRITE (6,*)  •  NOW  IN  HGPS2  SUBROUTINE' 
H ( 3 ) — 1.0D0 

MEAS.  NOISE 

R«DSQRT(SIGGPS**2+SIGCLK<2) **2 ) 

RETURN 


****************£ND  SUBROUTINE  HGPS2 ************************* 


END 


SUBROUTINE  HGPS3(H,R) 

THIS  SUBROUTINE  COMPUTES  THE  H  MATRIX  FOR  GPS  CLOCK  MEAS. 
GPS  -  CS017 

REAL* 8  H(10) , R, SIGGPS , SIGCLK ( 3 ) 

COMMON/MEASIG/SIGGPS , SIGCLK 
COMPUTE  THE  NON-ZERO  ENTRIES  OF  H 

WRITE (6,*)  '  NOW  IN  HGPS3  SUBROUTINE' 

H(5)— 1.0D0 

MEAS.  NOISE 

R=DSQRT (SIGGPS**2+SIGCLK (3) **2) 

RETURN 


*****************END  SUBROUTINE  HGPS3************************* 

END 

SUBROUTINE  HCLK1 (H, R) 

THIS  SUBROUTINE  COMPUTES  THE  H  MATRIX  FOR  ONLINE  CESIUM  CLOCK  MEAS. 

REAL* 8  H(10) ,R, SIGGPS, SIGCLK(3) 

COMMON/MEASIG/SIGGPS , SIGCLK 
COMPUTE  THE  NON-ZERO  ENTRIES  OF  H 
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WRITE (6,*)  '  HOW  IN  HCLK1  SUBROUTINE' 
H(l)  «  1„ 000 
H(3)  — 1.0D0 

MEAS.  NOISE 

R«DSQRT(SIGCLK(1) **2+SIGCLX(2) **2) 
RETURN 


END  SUBROUTINE  HCLK1 


END 

SUBROUTINE  HCLK2 (H, R) 

THIS  SUBROUTINE  COMPUTES  THE  H  MATRIX  FOR  ONLINE  CESIUM  CLOCK  MEAS. 

REAL*8  H(10) , R, SIGGPS , SIGCLK( 3 ) 

COMMON/MEASIG/SIGGPS , SIGCLK 
COMPUTE  THE  NON-ZERO  ENTRIES  OF  H 

WRITE (6,*)  '  NOW  IN  HCLK2  SUBROUTINE' 

H(l)  «  1.0D0 
H ( 5)  — 1.0D0 

MEAS.  NOISE 

R-DSQRT(SIGCLK(1) **2+SIGCLK(3 ) **2 ) 

RETURN 


*****************END  SUBROUTINE  HCLK2 ************************* 


END 

SUBROUTINE  KFUSER (MJDAY , TIME , DT , MODI , MOD2 , MOD3 , P , F , PHI , Q , QC , 

*  H,K,R,WT,W1,W2,W3,W4,PARMS,M,N,X) 

W3  ARRAY  CONTAINS  MEASUREMENT  RESIDUALS 
W4  ARRAY  CONTAINS  RAW  MEASUREMENTS 

REAL* 8  P(10,10) ,F( 10,10) , PHI (10,10) ,Q( 10,10 ) ,QC( 10,10 ) , 

*  H  ( 10) ,K( 10) ,WT( 10,10) ,W1( 10,10 )  ,W2 (10,10) ,X(10) ,SIG(10) 

*  , STIME, TIME, DT,R,R2D, MJDAY, W3 (10) ,W4 (10) 

INTEGER  PARMS (80), MODI , MOD2 , MOD3 

WRITE (6,*)  'NOW  IN  KFUSER,  PRINTING  P' 

WRITE (6 , *)  TIME 

WRITE ( 6 , * )  (P(I, I) ,1-1,10) 
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C 

DO  10  I  -  1 ,N 

SIG(I)  -  DSQRT ( P ( I , I ) ) 
10  CONTINUE 


WRITE {6,*)  ‘CLOCK  UNCERTAINTY  (1  SIGMA, micro-SEC) • 

WRITE ( 6 , * )  (SIG(I) ,1-1, N) 

STIME  -  TIME  -  0.5D0 
STIME  -  TIME 

IF (MODI .  EQ  •  *  LTUS ' )  CALL  LOTUS (X, SIG,W3 ,W4  ,MJDAY, STIME, N) 

RETURN 

END 

**************  end  OF  KFUSER  ********************* 


SUBROUTINE  LOTUS ( X , SIG , W3 , W4 , MJDAY , TIME , N) 


REAL*8  SIG (10), TIME 

REAL* 8  X(10) , MJDAY, W3 (10) ,W4(10) 

INTEGER  N 


WRITE(20, 100)  MJDAY, TIME, (X(J) ,J«1,N) , (SIG (I) , I«1,N) , 

+  (W3(K) ,K-1,5) , (W4(K) ,K«1,5) 

100  FORMAT ( IX ,  F7 . 1 ,  IX ,  F5 . 1 , 12  ( IX ,  F8 . 5 )  ,  5  ( IX ,  F8 . 5)  ,  5  ( IX ,  F8 . 5)  ) 

RETURN 

END 


SUBROUTINE  KFTABL ( EXTAB , I PTAB ) 

DIMENSION  ETAB (50,8) ,ITAB(50,8) , IPTAB ( 50 , 80) 
CHARACTER* 1  IREC(80) 

CHARACTER*4  I EVENT 
CHARACTER* 4  IPARMC 
CHARACTER* 4  IRESL 
DIMENSION  EXTAB (5 0,8) 

LOGICAL  EOF 
DIMENSION  I EVENT (8) 

EQUIVALENCE ( ETAB , ITAB ) 

I EVENT ( 1) 

I EVENT ( 2 ) 

I EVENT ( 3 ) 

I EVENT (4) 

IEVENT(5) 

IEVENT(6) 

IEVENT(7) 

I EVENT ( 8 ) 


i prop  1 
•UPDT' 
' PRNT* 
* STOR ' 
•RDUP' 
•USER' 
•STOP’ 
•TRAJ’ 


C 
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IPARMC  -  'PARM' 
EOF  «  .FALSE. 


ICNTE  -  0 
NREC  -  0 
ICNTP  -  0 
I BAD  -  0 
IPRCNT  ■  0 
IVNT  -  0 
HTRY  -  0 
MTRY  -  0 

WRITE  OUT  A  HEADER 


OPEN  THE  FILE  FOR  GETREC 
LUN-21 

OPEN  ( LUN ,  FILE- '  KFTABL  .IN',  ACCESS-  *  SEQUENTIAL  •  , 

1  FORM- • FORMATTED • ) 

WRITE(6, 90500) 

READ  A  RECORD 

10  CALL  GETREC ( LUN, IREC, NREC, EOF) 

CHECK  FOR  EOF 

IF (EOF)  GO  TO  9999 

CHECK  FIRST  4  CHARACTERS 

IFIRST  »  1 
ISIZE  -  4 

CALL  BLDWRD(IREC, IFIRST, ISIZE, IRESL) 

IS  IT  A  PARM  CARD? 

IF (IRESL  .EQ.  IPARMC)  GO  TO  6000 

NOT  A  PARM  SO  IS  IT  A  VALID  EVENT? 

ITYPE  -  0 
DO  100  1-1,8 

IF (IRESL  .EQ.  IEVENT(I) )  ITYPE  -  I 

100  CONTINUE 

WRITE (6, 101) ITYPE 

101  FORMAT  (IX,  ’ITYPE-  \I5) 

IF (ITYPE  .EQ.  0)  GO  TO  7000 

A  GOOD  EVENT  SO  PROCESS  IT 

FOR  NO  VALIDATION  OF  MODS  EXCEPT  FOR  SI/ZE  AND  NUMBER 
WE  WILL  USE  THIS  COMPUTED  GOTO  LATER 

GO  T0(5100, 5200, 5300, 5400, 5300, 5300, 5300, 5300) , ITYPE 
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GO  TO  5300 
C 

5300  CONTINUE 
C 

C  NO  LIMIT  ON  THE  PARMS  FIELDS  HERE 
C 

IF ( ITYPE  .EQ.  1)  GO  TO  5310 
NTRY  -  NTRY  +  1 
MTRY  -  NTRY 
ITAB(MTRY, 1)  -  IRESL 
C 

GO  TO  5320 

5310  IPRCNT  -  IPRCNT  +  1 
MTRY  -  IPRCNT  +  20 
ITAB(MTRY, 1)  -  IRESL 

5320  IF (NTRY  .GT.  20  .OR.  IPRCNT  .GT.  5)  GO  TO  7100 
C 

CALL  SGEN ( IREC , MTRY , ISTAT , ETAB) 

IF(ISTAT.EQ.l)  GO  TO  5330 
IVNT  -  1 
GO  TO  10 

5330  I BAD  -  I  BAD  +  1 
IVNT  -  0 
GO  TO  10 
C 

C  HERE  TO  PROCESS  PARMS 
C 

6000  CONTINUE 

IF (IVNT  .EQ.  0)  GO  TO  7200 
ICNTP  -  ICNTP  +  1 
IVNT  -  0 
DO  6100  1-1,80 

IPTAB ( ICNTP , I )  «  IREC (I) 

6100  CONTINUE 
C 

C  UPDATE  THE  POINTER  TO  THE  PARM  IN  THE  EVENT  RECORD 
C 

ETAB (MTRY, 8)  -  ICNTP 
C 

WRITE (6, 888) ICNTP 
888  FORMAT (IX, 'MAIN  IPTAB',15) 

C 

C 

GO  TO  10 
C 

C  OUTPUT  ERROR  MESSAGES 
C 

7000  WRITE (6, 90000)  IREC 
IBAD  *  IBAD  +  1 
IVNT  -  0 
GO  TO  10 

7100  WRITE(6, 90010) 

IBAD  -  IBAD  +  1 
IVNT  «  0 
GO  TO  10 

7200  WRITE (6, 90020)  IREC 
IBAD  -  IBAD  +  1 
IVNT  -  0 
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GO  TO  10 
C 

9999  CONTINUE 

IPCIBAO  .EQ.  0)  GO  TO  10000 
WRITE { 6 , 9003 0 ) IBAD 
STOP  77 

10000  CONTINUE 

WRITE(6, 90040) 

C  RETURN 

DO  3333  1-1,25 

WRITE (€ , 90100) ITAB(I , 1) , ITAB(I , 2) , ITAB(I , 3 ) , ITAB (I , 4 ) , 
1  ETAB (1,5) , ETAB (1,6) ,ETAB(I,7) ,ETAB(I,8) 

3333  CONTINUE 

IF  (ICNTP.GT.O)  THEN 
DO  3334  I— 1 , ICNTP 

WRITE (6, 90200) (IPTAB(I,J) ,J-1,80) 

3  334  CONTINUE 
END  IF 
C 

C  SWAP  THE  OUTPUT  TABLE  BACK 
C 

DO  10020  1-1,25 
DO  10010  J-1,8 

EXTAB ( I ,  J )  -  ETAB(I , J) 

10010  CONTINUE 
10020  CONTINUE 
C 

RETURN 

C 

C  FORMAT  STATEMENTS 
C 

90100  FORMAT (4 ( IX, A4) ,3(1X,F10.2) ,1X,F6.2) 

90200  FORMAT ( IX , 8 0A1 ) 

90000  FORMAT (IX, 'BAD  CARD  TYPE ' / IX , 8 0A1// ) 

90010  FORMAT (IX, 'TOO  MANY  PARM  OR  EVENT  CARDS'//) 

90020  FORMAT (IX, 'PARM  WITHOUT  EVENT'/IX, 80A1) 

90030  FORMAT (IX, 'UNSUCCESSFUL  TABLE  BUILD'/ 

1  IX, 'NUMBER  OF  ERRORS  DETECTED  -  ',110) 

90040  FORMAT (IX, 'SUCCESSFUL  TABLE  BUILD',//) 

90500  FORMAT (//,5X, 'EVENT  FILE',/) 

END 

SUBROUTINE  GETREC(LUN, INDAT,NREC, EOF) 

C 

CTHIS  ROUTINE  READS  IN  A  RECORD  FROM  (LUN)  AND  RETURNS 
C  THE  RESULT  AS  AN  80  CHARACTER  ARRAY,  INDAT 
C  IF  THE  READ  IS  SUCESSFUL, IRSTAT  -  1. 

C  IF  EOF  OCCURS,  IRSTAT-0 
C  IF  ANY  ERROR  OCCURS  IRSTAT— 1 
C  (NREC) -NUMBER  OF  RECORDS  READ 
C 

CHARACTER* 1  INDAT (80) 

LOGICAL  EOF 
C 

C  WRITE ( 6 , 888 ) 

888  FORMAT (IX, 'ENTER  GETREC') 

C 

c 
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C 

R£AD(UJN, 90000 , END-10) INDAT 
HREOHREC  +  1 
90010  FORXAT (IX, 80A1) 

EOF  ENCOUNTERED 

GO  TO  20 
10  EOF  -  .TRUE. 

NREC  -  NREC  -  1 

20  IF ( .NOT.  EOF)  WRITE(6, 90010) INDAT 
RETURN 

1000  FORMAT (80A1) 

END 

SUBROUTINE  FINDC ( ICH , I DATA , ISTART, IL0C) 

THIS  SUBROUTINE  SEARCHES  FOR  THE  CHARACTER  (ICH)  IN  THE 
80  CHARARCTER  STRING  (IDATA)  BEGINNING  AT  POSITION  (ISTART) 
AND  RETURNS  THE  LOCATION  OF  (ICH)  IN  (IL0C)  OR  RETURNS 
(ILOOO)  IF  (ICH)  IS  NOT  FOUND 

CHARACTER* 1  ICH, IDATA(80) 

MINOR  ERROR  CHECKING 

IF (ISTART  .LT.  1  .OR.  ISTART  ,GT.  80)  STOP  30 
IL0C  »  0 

DO  100  I  -  ISTART, 80,1 

IF(ICH  .NE.  IDATA (1) )  GO  TO  100 

FOUND  IT 

IL0C  -  I 
GO  TO  110 
100  CONTINUE 

110  CONTINUE 

WRITE ( 6 , 888 ) ICH , ISTART , ILOC 
888  FORMAT ( IX , 'FINDC' , IX, Al, 2 (IX, 13) ) 

RETURN 


END 

SUBROUTINE  BLDWRD ( IDATA , ISTART , ILEN , I OUT ) 

THIS  SUBROUTINE  TAKES  INDIVIDUAL  ARRAY  ELEMENTS  FROM 
(IDATA)  BEGINNING  AT  ( IDATA ( ISTART) )  FOR  LENGTH  (ILEN) 
CHARACTERS,  AND  ENCODES  THEM  INTO  A  SINGLE  WORD  IN  (IOUT) 
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C  A  MAXIMUM  OF  4  CHARACTERS  MAY  BE  ENCODED.  THIS  IS  FOR 
C  ALPHA  STRINGS  ONLY. 

C 

CHARACTER* 1  IDATA(80) 

CHARACTER* 4  I OUT 
IOUT-'  • 

C 

C  WRITE (6, 888) ISTART, ILEN 

888  FORMAT (IX, 'ENTER  BLDWRD  1ST  IEND' , IX, 2110) 

C 

C 

IEND  -  ISTART  +  ILEN  -  1 
GO  TO  (100,200,300,400) , ILEN 
C 

100  WRITE (IOUT, 90000)  I DATA (I END) 

GO  TO  10000 

200  WRITE (IOUT, 90010) (IDATA(I) ( I- I START , IEND) 

GO  TO  10000 

300  WRITE(IOUT, 90020) (IDATA(I) , I-ISTART, IEND) 

GO  TO  10000 

400  WRITE (IOUT, 90030) (IDATA(I) , I-ISTART, IEND) 

C 

C 

889  FORMAT ( IX ,  'BLDWRD  WORDOUT'  ,1X,  '  ('  ,A4,  ')  ') 

C 

C 

C10000  WRITE (6, 889) IOUT 
10000  CONTINUE 
RETURN 
C 

90000  FORMAT (Al) 

90010  FORMAT (2A1) 

90020  FORMAT (3A1) 

90030  FORMAT (4A1) 

END 

C 

SUBROUTINE  BLDNUM ( IDATA ,  ISTART,  IDEC,  IEND,  RESL) 

C 

C  ISTART  -  FIRST  DIGIT  OF  FIELD 
C  IEND  -  LAST  DIGIT  OF  THE  FIELD 
C  IDEC  «  LOCATION  OF  DECMAL  POINT 
C  IDATA  -  INPUT  STRING 
C  RESL  -  THE  RESULTING  NUMBER 
C 

C  CHARACTER* 10  ET 

CHARACTER*!  IDATA (80) 

C 

C  BUILD  INTEGER  PORTION  OF  THE  NUMBER 
C 

C  WRITE(6, 888) ISTART, IDEC, IEND 

888  FORMAT (IX, 'ENTER  BLDNUM' , 3 ( IX, 15) ) 

C 

C 

J  -  0 

RESL  -  0.0 
I LAST  -  IDEC  -  1 
ICNT  -  I LAST  -  ISTART 


C 


GPSSXM. FOR  APRIL  21,  1990  PAGE:  28 

DO  100  I  -  XSTART,ILAST, 1 

READ (IDATA (I) ,90000)  IDIGIT 
IPOW  -  10  **  (ICNT  -  J) 

RESL  -  RZSL  +  IDIGIT  *  IPOW 
J  -  J  +  1 
100  CONTINUE 
C 

C  NOW  DO  THE  FRACTIONAL  PORTION 
C 

C  IF  I END  -  IDEC  YOU  ARE  DONE 

C 

IF (IEND  .EQ.  IDEC)  GO  TO  1000 
C 

JSTART  -  IDEC  +  1 
K  -  1 

C 

DO  200  I  -  JSTART, IEND, 1 

READ(IOATA(I) ,90000)  IDIGIT 
C  WRITE (6,886)1, IDIGIT , IDATA ( I ) 

886  FORMAT ( IX , ' I , IDIG , IDAT  ',3110) 

XPOW  -  10.0  **  (-K) 

RESL  -  RESL  +  IDIGIT  *  XPOW 
K  •  K  +  1 
C 

C  WRITE(6, 887)1, IDIGIT, XPOW, RESL 

887  FORMAT ( IX , ' I-ID-IP-RESL* , 2 ( IX, 110)  ,  2 (IX, F14 . 8 ) ) 

C 

200  CONTINUE 
C 

1000  CONTINUE 
C 

C  WRITE (6, 889 )ISTART, IDEC, IEND, RESL 

889  FORMAT ( IX, 'BLDNUM  S  D  E  R' , 3 ( IX, 15) , F20 . 4 ) 


90000  FORMAT(Il) 

C 

END 

SUBROUTINE  SGEN ( IREC , NTRY , ISTAT , ETABLE ) 

DIMENSION  UtOD ( 3 ) 

DIMENSION  ETABLE (50, 8) 

CHARACTER* 1  IREC(80) 

CHARACTER*  1  ILPAR,  IRPAR,  ICOMM,  JA,  JB,  JF,  JT,  IBLK,  IPT 
CHARACTER* 4  AETAB, ISP 
ILPAR  -  ' ( ' 

IRPAR  -  ' ) ' 

ICOMM  -  ' , ' 

JA  -  'A' 

JB  «  'B' 

JF  -  'F' 

IBLK  »  '  ' 

IPT  -  ' • ' 

JT  ■  'T' 

ISP  -  '  ' 

C 

C  IREC-  INPUT  STRING 
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C  NTRY  -  THE  TABLE  ENTERY 
C  ISTAT-0  IP  OK,  -1  IF  ERROR 
C  LMOD  -  LOCATIONS  OF  COMMAS  IN  MODS 
C 

C  FIRST  SEARCH  FOR  MODS 
C  BEGIN  LOOKING  FOR  (  IN  CARD  COLUMN  5 
C 

C  WRITE ( 6 , 888) 

888  FORMAT (IX, 'ENTER  SGEN') 

C 

NMODS  -  0 
LMOD(l)  -  0 
IMOD(l)  -  0 
LMOD ( 2 )  -  0 
LMOD ( 3 )  -  0 
LMOD ( 2 )  -  0 
LMOD ( 3 )  -  0 
IREST  -  5 

CALL  FINDC(ILPAR, IREC, IREST, ILOC) 

IF  I  LOG  -  0  THEN  ASSUME  NO  MODS 
IF(ILOC  .EQ.  0)  GO  TO  1000 
FOUND  OPEN  LOOK  FOR  CLOSE 
IOP  -  ILOC  +  1 

CALL  FINDC(IRPAR, IREC, IOP, ILOC) 

IF  WE  DID  NOT  FIND  )  THEN  AN  ERROR 
IF (ILOC  .EQ.  0)  GO  TO  5000 

IOP  WILL  BE  THE  FIRST  CHARACTER  OF  MODS,  IOC  THE  LAST  CHARACTER 
NOT  INCLUDING  THE  PARENS 

IOC  -  ILOC  -  1 

NOW  SEE  HOW  MANY  COMMAS  IN  THE  () 

BETTER  NOT  BE  MORE  THAN  2 

NMODS  -  1 

DO  10  I  -  IOP, IOC, 1 

IF(IREC(I)  .NE.  ICOMM)  GO  TO  10 
LMOD (NMODS)  -  I 
NMODS  -  NMODS  +  1 
IF (NMODS  .GT.  3)  GO  TO  5000 
10  CONTINUE 
C 

C  NOW  IOP  IS  THE  FIRST  CHARACTER  IN  THE  STRING  OF  MODS 

•C  IOC  IS  THE  LAST  CHARACTER 

C  LMOD ( 1 ) ( 2 )  CONTAINS  EITHER  0  OR  THE  LOCATION  OF  THE  COMMAS 
C  NMODS  IS  THE  NUMBER  OF  MODS 

C  FOR  NOW  MODS  ARE  PROCESSED  DUMB  WILL  FIX  LATER 
C 

C  DO  MOD  1 

IREST  -  IOC  +  1 
IFIRST  -  IOP 
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I LAST  -  IOC 

ILEN  -  IMOD ( 1 )  -  I FIRST 

IF (NMOOS  .EQ.  1)  ILEN  -  I LAST  -  IFIRST  +  1 
CALL  BLOMRD ( IREC , IFIRST , ILEN , AETAB) 

STABLE (NTRY, 2)  -  AETAB 
C 

C  IS  THERE  ANOTHER  MOD? 

C 

IF (NMODS  .EQ.  1)  GO  TO  1000 
C 

C  DO  MOD  2 
C 

ILEN  -  IMOD ( 2 )  -  LMOD(l)  >1 

IF (NMODS  .EQ.  2)  ILEN  -  IOC  -  LMOD(l) 

IFIRST  -  IMOD(l)  +  1 

CALL  BLDWRD(IREC, IFIRST, ILEN, AETAB) 

STABLE  (NTRY,  3)  -  AETAB 
C 

IF (NMODS  .EQ.  2)  GO  TO  1000 
C 

C  DO  THE  THIRD  MOD 
C 

ILEN  -  I LAST  -  LMOD ( 2 ) 

IFIRST  -  LMOD ( 2 )  +  1 

CALL  BLDWRD ( IREC , IFIRST , ILEN , AETAB ) 

ET ABLE (NTRY , 4 )  -  AETAB 

C 

C  MODS  ARE  DONE 
C 

1000  CONTINUE 
C 

C  MODS  ARE  DONE,  NOW  DO  THE  REST 

C  (IREST)  CONTAINS  THE  STARTING  LOCATION  FOR  THE  REST  OF  THE  CARD 
C  FIRST  SEE  IF  IT  IS  AN  (AT)  OR  A  (FROM) 

C 

CALL  FINDC(JA, IREC, IREST, ILOC) 

IF(ILOC  .EQ.  0)  GO  TO  2000 
C 

C  PROCESS  THE  (AT)  TIME 
C 

ETABLE  ( NTRY , 6 )  -  0 
ET  ABLE (NTRY, 7)  »  0 
C 

C  ILOC  CONTAINS  THE  POSITION  OF  THE  (A) 

C  ILOC  +  3  SHOULD  BE  THE  FIRST  DIGIT 
C 

11  -  ILOC  +  3 
C 

C  FIND  THE  END  POINT  OF  THE  NUMBER 
C 

CALL  FINDC(IBLK, IREC, II, ILOC) 

IF ( ILOC  .EQ.  0)  GO  TO  5000 

12  «  ILOC  -1 
C 

C  NOW  FIND  THE  LOCATION  OF  THE  (DECIMAL  POINT) 

C  PUT  IT  IN  ( IPD) 

C 

CALL  FINDC (IPT, IREC , II , IPD) 
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IF ( IPD  .EQ.  0  .OR.  IPD  .GT.  ILOC)  GO  TO  5000 
C 

C  IPD  -  LOCATION  OF  DECIMAL 
C  II  -  FIRST  DIGIT 
C  12  -  LAST  DIGIT 
C  BUILD  ETAB(NTRY ,5} 

C 

CALL  BLDNUM(IR£C, II, IPD, 12 , ETABLE (NTRY, 5) ) 

THE  (AT)  CARD  IS  DONE,  SET  THE  FLAG  GOOD  AND  RETURN 

ISTAT  -  0 
GO  TO  9000 

2000  CONTINUE 

THIS  IS  A  (FROM)  (TO)  (BY)  CARD 

FIND  WORD  (FROM) 

CALL  FINDC(JF, IREC, IREST, ILOC) 

IF  (ILOC  .EQ.  0)  GO  TO  5000 
ASSUME  BLANKS  ON  BOTH  SIDES  OF  NUMBERS 

11  «  ILOC  +  5 

HOW  MANY  DIGITS? 

CALL  FINDC(IBLK, IREC, II, ILOC) 

IF (ILOC  .EQ.  0)  GO  TO  5000 

12  -  ILOC  -  1 

FIND  THE  DECIMAL  POINT 

CALL  FINDC(IPT, IREC, II, ILOC) 

IF (ILOC  .EQ.  0  .OR.  ILOC  .GT.  12)  GO  TO  5000 
IPD  -  ILOC 

CALL  BLDNUM (IREC, II, IPD, I 2, ETABLE (NTRY, 5) ) 

NOW  DO  THE  (TO)  WORD 
ISTART  -  12 

CALL  FINDC(JT, IREC, ISTART, ILOC) 

IF (ILOC  .EQ.  0)  GO  TO  5000 

SHOULD  BE  A  BLANK  BETWEEN  (TO)  AND  NUMBER 

11  -  ILOC  +  3 

FIND  12,  THE  LAST  DIGIT 

CALL  FINDC(IBLK, IREC, II, ILOC) 

IF (ILOC  .EQ.  0)  GO  TO  5000 

12  «  ILOC  -  1 

FIND  THE  PERIOD 
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CALL  FINDC ( IPT , IREC , II, ILOC) 

IF(IL0C  .EQ.  0  .OR.  ILOC  .GT.  12)  GO  TO  5000 
IPO  -  ILOC 

CALL  BL0NUM(IREC,X1,XPD,X2,ETABLE(NTRY,6) ) 

DO  THE  (BY) - FIND  B 

ISTART  -  12 

CALL  FINDC(JB/ IREC, ISTART, ILOC) 

IF (ILOC  .EQ.  0)  GO  TO  5000 

IF  YOU  GET  HERE  YOU  GOT  (B)  SO  GET  A  NUMBER 

11  -  ILOC  +  3 

FIND  THE  LAST  DIGIT 

CALL  FINDC(IBLK, IREC, II, ILOC) 

IF (ILOC  .EQ.  0)  ILOC  -  80 

12  -  ILOC  -  1 


FIND  THE  DECIMAL  POINT 

CALL  FINDC (IPT, IREC, II, ILOC) 
IF (ILOC  .EQ.  0  .OR.  ILOC  .GT. 
IPD  -  ILOC 


12)  GO  TO  5000 


DO  IT 

CALL  BLDNUM(IREC, II, IPD, 12 , ETABLE (NTRY, 7) ) 

A  DEEP  SIGHKHHHHHHH,  WE  ARE  DONE! !!!!!! 

ISTAT  -  0 
GO  TO  9000 


WE  HAVE  A  PROBLEM  SO  COME  HERE 


5000  CONTINUE 

WRITE (6, 90000)  IREC 
ISTAT  -  1 
GO  TO  9000 


9000  RETURN 


FORMAT  STATEMENTS 


90000  FORMAT (IX, ’ERROR  IN  FOLLOWING  CARD'/80A1//) 
END 
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INTRODUCTION 


This  User’s  Guide  provides  instructions  for  using  a  software  package  hereafter  referred 
to  as  the  Covariance  Simulation  Program  (CSP).  This  guide  is  intended  to  familiarize  the  user  with 
the  CSP  operation.  The  hardware  suite  needed  to  operate  the  program  is  listed  and  detailed 
instructions  for  its  use  are  provided. 

The  CSP  was  developed  by  SYNETICS  for  the  U.S.  Coast  Guard  Omega  Navigation 
System  Center  under  contract  DTCG-23-86-A-20022.  The  models  used  in  this  simulation  are 
described  in  a  report  entitled  "Analysis  of  GPS  Timing  Data  in  Support  of  Omega  System 
Synchronization:  A  Cesium  Stability  Study",  also  developed  under  the  same  contract.  The  report 
provides  details  of  the  simulation,  a  source  code  listing  (addendum),  and  includes  a  discussion  of 
simulations  based  on  timing  data  from  Omega  Stations  Norway  and  Hawaii. 


USING  THE  PROGRAM 


This  section  of  the  report  describes  the  hardware  requirements,  installation  process,  and 
operational  procedures  required  to  use  the  CSP.  These  instructions  will  allow  a  user  who  is  familiar 
with  the  format  and  contents  of  the  input  data  to  successfully  operate  the  program. 

2.1  HARDWARE  CONFIGURATION  ASSUMPTIONS 

The  CSP  is  distributed  on  one  5.25  inch  high  density  diskette.  This  diskette  contains 
an  automatic  installation  procedure  which  will  load  all  necessary  files  to  the  user’s  hard  disk  drive. 

The  user  computer  is  assumed  to  be  configured  as  follows: 

•  IBM  PC- AT  or  AT-compatible  clone 

•  Hard  Disk  Drive,  configured  as  Logical  Device  C:  with  at  least  1  megabyte  of 
free  RAM  memory 

•  80287  Math  Co-Processor 

•  1.2  megabyte  5.25-inch  floppy  disk  drive 

•  Version  3.1  (or  higher)  DOS  or  PC-compatible  operating  system,  such  as 
Microsoft  DOS  or  XENIX,  resident  on  the  hard  disk 

•  LOTUS  1-2-3  Release  2.01  installed  in  a  subdirectory  named 
C:\LOTUS  on  drive  C: 

•  A  color  or  monochrome  graphics  adapter  and  display. 

The  CSP  package  will  automatically  install  and  can  be  operated  on  a  system  which  meets  these 
minimum  criteria. 
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INSTALLATION 


The  following  steps  will  allow  you  to  automatically  install  all  necessary  program  and 
data  files  for  operation  of  the  CSP. 


1.  Turn  on  your  computer. 

2.  Once  you  are  at  your  DOS  prompt  (ex.  Cj\)  insert  the  program  diskette  in  dri^e  A. 

3.  Type  the  string  A:  followed  by  ENTER.  You  should  now  see  the  DOS  prompt 
on  your  screen. 

4.  Enter  the  command  INSTALL  followed  by  ENTER. 


The  installation  procedure  will  create  a  subdirectory  named  C:\FILTER  and  all  necessary  program 
and  data  files  will  be  copied  to  that  directory.  While  transferring  the  files  from  the  program 
diskette  to  the  subdirectory,  the  installation  procedure  lists  the  file  names  on  the  screen  for  the 
user’s  information.  At  the  end  of  the  installation  process,  the  following  message  will  appear  to 
indicate  successful  installation  and  provide  the  necessary  user  instructions  to  begin  the  execution 
of  the  system: 


INSTALLATION  COMPLETE 
TYPE  C: 

THEN  CD  \FILTER 
THEN  FILTER 

To  begin  execution  of  the  system 


Appendix  A  lists  the  directories  (and  their  contents)  resulting  from  the  installation 
procedure.  Once  the  installation  is  complete,  you  are  ready  to  begin  using  the  system. 
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SYSTEM  SETUP 


You  can  initiate  the  execution  of  the  CSP  at  any  time  by  entering  the  following 
commands  from  the  DOS  prompt. 

1.  Type  Cx  followed  by  an  ENTER. 
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2. 

3. 


At  the  C:\  prompt  type  CD  VFILTER  followed  by  an  ENTER. 
Type  FIL  TER  followed  by  an  ENTER. 


A  Filter  Menu  (Figure  2.3-1)  will  appear  on  the  screen.* 


FILTER  MENU 


1  -  TRANSLATE  FILE 

2  -  EDIT  EVENTS  FILE 

3  -  EDIT  PARAMETER  FILE 

4  -  RUN  A  FILTER 

5  -  LOTUS 

6  -  QUIT 


January  31,  1990  10:33:57  am 


CAPS 


Memory:  488  K 


Press  H  for  Hslp 


Figure  23-1  Filter  Menu  Selection  Screen 

The  menu  selection  is  made  in  one  of  several  wavs: 

•< 

1.  Move  the  arrow  pointer  to  point  to  your  selection .  and  press  ENTER. 

You  move  the  arrow  pointer  by: 

■  pressing  the  space  bar  to  move  it 
downward  through  the  menu  selection 

*  Filter  Menu  was  generated  using  the  Automenu  Software  Management  System.  This  is  provided 
with  the  CSP  package. 
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-  pressing  the  up  or  down  arrow  to  move  the 
pointer  up  or  down  through  the  menu 
selections 

2.  Highlight  your  menu  selection  by  the  same  method  as  shown  above. 

3.  Make  your  selection  directly  by: 

Pressing  the  number  key  or  function  key  that  corresponds  to  your  selection 
number 


The  Filter  Menu  is  accompanied  by  the  on-line  help  feature  that  is  activated  by  typing 
H  for  the  "menu  help  key"  display.  Typing  H  again  returns  an  extended  on-line  help  narrative. 
Execution  of  the  CSP  Allows  steps  1-6  in  the  Filter  Menu.  Once  the  execution  of  the  previous 
selection  is  complete,  the  user  is  brought  back  to  the  Filter  Menu  selection  screen.  The  following 
sections  discuss  the  operation  of  each  step  in  detail. 


2.4  FILE  TRANSLATION 

The  first  step  in  running  the  system  is  the  translation  process.  This  is  Menu  Option  1 
(Figure  2.3-1).  The  translation  process  prepares  USNO  Remote  Data  Acquisition  System  (RDAS) 
data  (see  Section  3.2)  for  input  into  the  filter.  This  step  has  to  be  executed  first,  any  time  a  new 
set  of  timing  data  is  collected  and  the  user  wishes  to  use  this  data  as  an  input  to  CSP.  Copy  the 
file  of  USNO  RDAS  data  into  the  filter  subdirectory  using  the  DOS  COPY  command.  Then  use 
program  FILEPROC,  activated  by  the  above  Filter  Menu  selection,  to  scan  the  input  data  file  and 
retrieve  the  required  six  record  fields  as  indicated  in  Tables  3.2-1  and  3.2-2,  respectively.  These 
data  must  be  configured  as  shown  in  Figure  3.2-1.  At  this  point,  it  should  be  noted  that  the 
program  diskette  contains  timing  data  files  NORWAYl.DAT  and  HAWAII1.DAT  for  Omega 
Stations  Norway  and  Hawaii,  respectively. 

When  Option  1  is  selected  from  the  menu,  the  following  messages  will  be  displayed  on 

the  screen: 

TRANSLATING  FILE 

PLEASE  ENTER  INPUT  USNO  RAW  DATA  FILENAME: 
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In  response  to  the  filename  question,  enter  the  name  assigned  to  the  USNO  PTTI  data 
file  which  was  copied  into  the  FILTER  subdirectory  for  processing.  The  input  data  file  to  be 
translated  can  have  anv  valid  DOS  file  name.  For  this  example  we  are  using  the  file  named 
NORWAYl.DAT  provided  with  CSP  package. 

The  translation  process  creates  an  output  file  specified  by  the  user  at  the  "PLEASE 
ENTER  OUTPUT  DATA  FILENAME:"  prompt.  This  file  will  be  used  as  input  for  the  filter 
execution  which  is  activated  by  Menu  Option  4.  The  format  of  the  translation  output  file  is  shown 
in  Table  3.2-3  for  Norway  timing  data  and  in  Table  3.2-4  for  Hawaii  timing  data.  For  this  example 
NORWAY.DAT  filename  is  used. 

Within  the  translation  process  is  a  conversion  of  RDAS  timing  data.  The  RDAS  data 
is  in  units  of  seconds  and  in  the  format  described  in  Section  3.3  of  reference  1.  For  use  in  the  CSP 
it  is  necessary  to  "convert"  the  data  by  scaling  it  to  microseconds  and  changing  the  format.  This 
conversion  is  accomplished  by  a  simple  algorithm  described  below: 

1.  a  raw  measurement  is  read  from  an  input  data  file 

2.  if  a  measurement  is  Jess  than  0.5  seconds,  it  is  multiplied  by  106; 
otherwise,  1  is  subtracted  from  it  and  then  it  is  multiplied  by  106 

3. ’  the  converted  measurement  is  written  to  an  output  data  file 

The  translation  process  will  take  a  few  moments  to  complete.  At  termination,  the 
following  message  will  be  displayed: 

STOPPED  PROCESSING  NORWAY.DAT  FILE 
Stop  -  Program  terminated 


*  Although  the  raw  measurements  have  been  processed  as  described  they  will  be  referred  to  as  raw- 
measurements 


An  example  of  a  complete  translation  session  is  shown  below: 


TRANSLATING  FILE... 

PLEASE  ENTER  INPUT  USNO  (RDAS)  RAW  DATA  FILENAME  : 
NORWAYl.DAT 

PLEASE  ENTER  OUTPUT  DATA  FILENAME  : 

NORWAY.DAT 

STOPPED  PROCESSING  NORWAY.DAT  FILE 
Stop  -  Program  terminated. 


At  the  end  of  the  translation  process,  the  user  is  automatically  brought  back  to  the  Filter  Menu 
selection  screen. 


2.5  EVENTS  FILE  EDITING 


The  diskette  contains  a  sample  events  file.  You  may  view  this  file  by  printing  file 
KFTABL.IN.  An  example  of  an  events  file  is  shown  in  Figure  2.5-1.  The  event  file  controls  the 
sequence  of  execution  of  the  CSP.  There  are  four  types  of  events  which  are  described  in  detail 
in  Section  3.1  (see  also  Reference  1).  Menu  Option  2  provides  the  user  with  the  capability  to 
modify  and  update  the  events  file  for  particular  applications. 


At  this  point,  the  user  has  to  make  a  decision:  to  customize  the  input  file  for  the 
specific  simulation  application  or  skip  the  editing  session  and  use  the  default  values.  Shareware 
program  PC- WRITE  is  provided  and  invoked  to  edit  the  events  file.  The  name  of  this  file  is  always 
KFTABL.IN.  When  first  invoked,  PC-WRITE  asks  if  you  wish  to  make  a  backup  copy  of 
KFTABL.IN.  If  you  do,  then  press  function  key  F9  If  not,  then  press  ESC.  At  any  time  during 
the  editing  process  you  can  invoke  on-screen  help  by  pressing  function  key  FI  twice. 
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PROP  (CNST,  STAT)  FROM  0.0  TO  31.0  BY  0.5 
UPDT  (VARY,  GPS1,  STAT)  FROM  0.0  TO  31.0  BY  3.0 
UPDT  (VARY,  GPS2,  STAT)  FROM  0.0  TO  31.0  BY  3.0 
UPDT  (VARY,  GPS3,  STAT)  FROM  0.0  TO  31.0  BY  3.0 
USER  (LTUS)  FROM  0.0  TO  31.0  BY  0.5 
STOP  AT  31.0 


Figure  2.5-1  Sample  Event  File 


Once  the  editing  process  is  complete,  save  the  KFTABL.IN  events  file  and  exit  PC- 
WRITE  by  pressing  function  keys  FI  then  F2  consecutively.  Please  note  that  you  do  not  need  to 
use  PC- WRITE  to  edit  the  events  file.  You  can  use  any  text  editor  which  creates  ASCII  files. 

When  the  events  file  editing  session  is  terminated,  the  user  is  automatically  returned 
to  the  Filter  Menu  selection  window. 


PARAMETER  FILE  EDITING 


Your  diskette  contains  a  sample  parameter  file.  You  may  view  this  file  by  printing  file 
P ARM. IN.  An  example  of  a  PARM.IN  file  is  shown  in  Figure  2.6-1.  The  values  in  the  parameter 
list  are  constants  that  describe  the  characteristics  of  the  covariance  "scenario”  to  be  simulated. 
Table  2.6-1  lists  the  parameter  names  and  units. 


The  clock  record,  designated  by  "CLK="  name,  contains  the  initial  phase  offset  and 
frequency  offset  (Ref.  1)  estimates  for  each  cesium  clock:  on-line,  secondary  and  primary  in  that 
order.  The  next  four  records  define  the  cesium  clock  error  model  and  the  corresponding  process 
noise  (Ref.  1).  The  final  record,  designated  by  "SIG=-"  name,  defines  the  uncertainties  for  GPS  and 
cesium  clock  time  measurements.  Menu  Option  3  provides  the  user  with  the  capability  to  modify 
and  update  the  parameter  file  as  desired.  At  this  point,  the  user  has  two  choices:  to  customize  the 
input  file  for  the  specific  simulation  application  or  skip  the  editing  session  and  use  the  default 
values.  Shareware  program  PC- WRITE  is  provided  and  invoked  to  edit  the  parameter  file.  The 
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’CLK 

TIME 

’COV 

’ERR 

QCOV 

’SIG 

I 

I 


Figure  2.6-1  Sample  Parameter  File  for  Norway  Data  Processing 


Table  2.6-1 

PARAMETER  FILE  CONTENTS 


=  ’,0.3570,  0.011461,  -0.7380,  0.119371,  1.21,  -0.050453 
=  ’,0.5,  0.5,  0.5 

=  ’,0.0001,  0.0016,  0.000225,  0.0016,  0.0000,  0.0016 
=  ’,0.03,  15.0,  0.03,  15.0,  0.03,  15.0 
=  ’,0.00036.  0.0,  0.00036,  0.0,  0.00036,  0.0 
=  ’,0.050,  0.007,  0.010,  0.008 


RECORD 

NAME 

VALUE  (UNITS) 

1 

CLK 

Initial  time  offset  estimate  (Atsec),  initial  frequency 
estimate  (/xsec/day). 

2 

TIME 

Paramter  used  in  calculation  of  Q  matrix  elements  (days) 

3 

COV 

Mean  square  uncertainty  in  phase  offset  estimate  (/u.sec)2, 
mean  square  uncertainty  in  frequency  offset  estimate 
(/xsec/day)2. 

4 

ERR 

Largest  jump  in  cesium  frequency  offset 

(/isec/day),  corresponding  to  measurement  interval  (days) 

5 

O  COV 

Mean  square  frequency  fluctuation  (^isec)2.  mean  square 
frequency-rate  fluctuation  (/xsec/day)2. 

6 

SIG 

Root  mean  square  time  measurement  error  (yusec). 

GPS  or  cesium  clocks 
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name  of  this  file  is  always  P ARM. IN.  When  first  invoked,  PC- WRITE  asks  if  you  wish  to  make  a 
backup  copy  of  PARM.IN.  If  you  do,  then  press  function  key  F9;  if  not,  then  press  ESC.  At  any 
time  during  the  editing  process  you  can  invoke  on-screen  help  by  pressing  function  key  FI  twice. 


Once  you  have  completed  the  editing  process,  save  the  PARM.IN  parameter  file  and 
exit  PC-WRITE  by  pressing  function  keys  FI  then  F2  consecutively.  Please  note  that  you  do  not 
have  to  use  PC-WRITE  to  edit  the  parameter  file.  You  can  use  any  text  editor  which  creates 
ASCII  files. 


2.7  FILTER  EXECUTION 

Once  you  have  created  your  input  files,  you  are  now  ready  to  run  the  filter.  The  filter 
will  accept  as  input  the  following  files: 

NORWAY.DAT  (or  HAWAII.DAT)  -  created  by  Menu  Option  1 
(translation  process) 

KFTABL.IN  -  created  by  Menu  Option  2 
PARM.IN  -  created  by  Menu  Option  3. 

The  above  files  are  automatically  opened  and  read  by  the  CSP.  As  output,  the  Kalman  filter 
portion  of  the  CSP  (GPSSIM)  creates  a  file  which  will  be  input  to  LOTUS  1-2-3.  LOTUS  is  used 
to  create  the  graphs  of  data  processed  during  early  steps  for  analysis  and  display  to  the  user. 

To  invoke  the  filter  process,  select  Option  4  from  the  menu.  The  following  will  be 

displayed: 

RUNNING  FILTER  PROGRAM... 

PLEASE  ENTER  RAW  MEASUREMENT  DATA  FILE 


In  response  to  the  question,  enter  the  name  you  have  selected  for  the  output  file  of  the  translation 
process.  Next,  the  user  has  to  enter  the  name  of  a  filter  output  file:  for  example,  you  might  select 
LOTUS.OUT  as  an  output  file  name.  Figure  2.7-1  shows  an  example  of  a  LOTUS  spreadsheet 
generated  by  importing  this  file.  Table  2.8-1  defines  the  data  elements  in  the  output  file.  A  file 
import  and  other  LOTUS  functions  are  described  in  some  detail  in  the  next  section.  The  filter 


output  file  is  an  ASCII  file  and  can  be  listed  for  review  prior  to  the  next  step  if  desired. 
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The  filter  execution  process  will  take  a  minute  or  two  to  complete  depending  on  the 
entries  such  as  the  start  and  stop  times  in  the  event  file.  At  termination,  the  following  messages 
will  be  displayed: 

Stop  -  Program  terminated. 

DIRECTORY  LISTING  OF  .OUT  FILES 

Volume  in  drive  C  has  no  label 
Directory  of  C:\FILTER 

LOTUS  OUT  size  mm-dd-yy  hh:mm 
Strike  a  key  when  ready... 

where  mm-dd-yy  and  hh:mm  fields  contain  the  date  and  time  of  file  generation. 
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Figure  2.7-1  Sample  Worksheet  Generated  from  LOTUS.OUT  File 
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Figure  2.7-1  (con’t)  Elements  I  through  P  of  the  Sample  Worksheet 
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Figure  2.7-1  (con’t)  Elements  Q  through  S  of  the  Sample  Worksheet 
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USING  LOTUS 


Menu  Option  5  will  invoke  LOTUS  1-2-3.  LOTUS  is  used  to  provide  a  data 
manipulation  and  graphing  capability  for  the  output  data  file.  The  following  discussions  will  assume 
that  you  have  created  a  file  named  LOTUS.OUT  as  the  output  of  the  filter  process. 


Select  Menu  Option  5  to  invoke  LOTUS  1-2-3.  At  the  main  LOTUS  prompt,  select 
1-2-3.  When  the  main  spreadsheet  screen  appears,  follow  the  steps  given  below  to  load  the  filter 
output  file  (LOTUS.OUT)  into  the  spreadsheet  for  processing.  The  following  commands  are 
provided  to  facilitate  use  of  LOTUS: 


1.  /  -  Enter  Lotus  command  mode 

2.  F  -  File  access  sub-menu 

3.1  -  Import  file 

4.  N  -  Import  file  as  numbers 


You  will  now  be  prompted  for  a  file  name.  Assuming  a  file  named  LOTUS.OUT,  press  ESC,  then 
type  the  file  name  C:\F1LTER\LQTUS.QUT  followed  by  the  ENTER  key  followed  by  the  ENTER 
key  again. 
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The  LOTUS.OUT  data  file  will  now  be  imported  into  your  spreadsheet  for  processing. 
This  file  contains  19  data  elements  for  each  time  point  as  shown  in  Table  2.8-1. 

Once  the  data  are  loaded  into  the  spreadsheet,  the  user  can  select  the  combination  of 
data  elements  to  be  displayed  graphically.  Only  data  elements  C  through  S  can  be  graphed.  The 
graphs  are  presented  in  standard  rectangular  xy-coordinates  with  time  on  the  positive  X  axis.  Any 
graphics  created  during  these  analyses  can  be  printed  using  the  PrintGraph  option  from  the  main 
LOTUS  1-2-3  menu.  The  user  performs  the  selection  in  the  following  manner: 

1.  Make  sure  the  main  LOTUS  1-2-3  menu  is  on  your  screen. 

2.  Move  the  menu  pointer  to  PrintGraph  and  press  ENTER. 

3.  Follow  the  instructions  on  the  screen  if  the  program  asks  you  to 
change  a  disk 

4.  Enter  the  command: 

I  -  Select  Image-Select.  Follow  the  instructions  on  the 
screen  to  indicate  the  graphs  you  want  to  print. 

A  -  Select  Align  to  tell  PrintGraph  that  the  paper  is  currently 
positioned  at  the  top  of  the  page. 

G  -  Select  GO  to  begin  printing. 

P  -  Select  Page  to  advance  the  page  when  printing  is 
complete. 

5.  Type  e  to  select  EXIT. 

6.  Type  y  to  confirm  that  you  want  to  leave  PrintGraph. 
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Table  2.8-1 

LOTUS  Output  Data  File  Contents 
(NORWAY) 


ELEMENT 

CONTENT 

A 

Data  Point  Modified  Julian  Date 

B 

Data  Point  Time 

C 

Phase  Offset  -  Online  Cesium  Clock 

D 

Frequency  Offset  -  Online  Cesium  Clock 

E 

Phase  Offset  -  Primary  Cesium 

F 

Frequency  Offset  -  Primary  Cesium 

G 

Phase  Offset  -  Secondary  Cesium 

H 

Frequency  Offset  -  Secondary  Cesium 

I 

Standard  Deviation  -  Element  C  above 

J 

"  "  *  Element  D  above 

K 

-  Element  E  above 

L 

-  Element  F  above 

M 

"  Element  G  above 

N 

-  Element  H  above 

O 

GPS  -  Online  Cesium  Measurement  Residual 

P 

GPS  -  Secondary  Cesium  Measurement  Residual 

Q 

GPS  -  Primary  Cesium  Measurement  Residual 

R 

Online  -  Secondary  Cesium  Measurement  Residual 

S 

Online  -  Primary  Cesium  Measurement  Residual 
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Once  you  have  completed  work,  the  spreadsheet  can  be  saved  to  disk  for  future  reference 
by  entering  the  following  commands  then  supplying  a  file  name: 

1.  /  -  Enter  LOTUS  command  mode 

2.  F  -  File  access  sub -menu 

3.  S  -  Save  file 

The  user  will  be  automatically  returned  to  the  Filter  Menu  after  exiting  LOTUS. 


COVARIANCE  SIMULATION 


3.1  GENERAL  INFORMATION 

The  CSP  is  used  to  predict  the  phase  and  frequency  offsets,  and  the  prediction 
uncertainties  for  the  three-cesium  frequency  standard  ensemble  at  an  Omega  transmitting  station. 
A  Kalman  filter  is  employed  to  propagate  a  phase  offset  and  a  phase  drift  rate  (frequency  offset) 
covariance  matrix,  and  to  update  the  estimation  error  covariance,  based  on  assumed  cesium 
standard  parameters  and  initial  offset  data.  Five  sets  of  timing  measurements  are  employed  in  the 
estimation  process: 

•  GPS  minus  on-line  cesium  clock 

•  GPS  minus  primary  cesium  clock 

•  GPS  minus  secondary  cesium  clock 

•  On-line  cesium  clock  minus  primary  cesium  clock 

•  On-line  cesium  clock  minus  secondary  cesium  clock 

There  are  three  types  of  outputs  generated  for  each  of  cesiums: 

•  estimates  of  phase  and  frequency  offsets 

•  filter  determined  uncertainties  in  these  estimates 

•  phase  offset  residuals  (difference  between  phase  offset  measurement  and  phase 
offset  estimate). 

Figure  3.1-1  illustrates  the  computational  flow  of  the  simulation.  The  values  in  the 
parameter  file  are  constants  that  describe  the  characteristics  of  the  simulated  system.  Table  2.6-1 
lists  the  parameter  names  and  units.  Each  line  in  the  table  is  a  record  in  the  P ARM. IN  file  (an 
example  of  a  PARM.IN  file  is  shown  in  Figure  2.6-1). 


3-1 


An  event  file  (KFTABLIN)  controls  the  sequence  of  execution  of  the  covariance 
analysis  software.  There  are  four  types  of  events  used  by  the  CSP: 


•  The  update  event  (UPDT)  updates  the  covariance  to 
incorporate  timing  measurements 

•  The  propagate  event  (PROP)  propagates  the  covariance 
through  a  time  interval 

•  The  user  event  (USER)  writes  data  to  the  filter  output  file 

•  The  stop  event  (STOP)  terminates  the  simulation. 


Events  are  scheduled  by  order  of  time  sequence  labeling  and  by  order  of  appearance 
in  the  event  list.  If  several  events  are  scheduled  at  one  time,  they  will  be  executed  in  the  order  they 
appear  in  the  event  table.  The  sole  exception  is  the  propagation  event  from  time  tk  to  tk+1  which 
is  defined  at  tk.  It  is  executed  after  all  other  events  at  tk  and  before  any  other  event  at  tk+ 1  . 

The  USER  event  is  a  general  purpose  exit  from  the  covariance  simulation  to  perform 
auxiliary  computations  on  the  covariance  analysis  in  progress.  In  the  simulation  it  is  used  to 
generate  a  LOTUS  output  file.  The  LOTUS  file  is  a  list  of  records,  each  one  containing  time, 
phase  and  frequency  estimates,  RMS  uncertainties  in  the  estimates  and  measurement  residuals. 

The  STOP  event  provides  the  program  the  ability  to  terminate  normally  without 
processing  an  entire  event  file.  No  records  after  the  first  STOP  are  processed.  This  permits  the 
user  to  modify  or  debug  a  simulation  easily.  A  full  discussion  of  the  simulation  model  is  contained 
in  Reference  1. 

The  list  below  describes  the  operational  flow  of  the  simulation  in  the  order  that 
corresponds  to  the  CSP  execution. 
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PHASE  AND 
FREQUENCY 
OFFSETS 


1.  Parameter  file  (PARM.IN)  is  read  containing  initial  phase  and 
frequency  offsets,  parameters  defining  the  cesium  clock  error 
model,  mean  square  uncertainties  in  phase  and  frequency  offset 
estimates  and  measurement  errors.  A  description  of  parameter 
file  contents  is  provided  in  Table  2.6-1. 

2.  The  constants  contained  in  the  parameter  file  are  used  to  set 
the  following  Kalman  filter  matrices: 

•  F  (transition  matrix)  -  transition  matrix  allows 
calculation  of  the  state  vector  (phase  and  frequency 
offset)  at  some  time  t,  given  complete  knowledge  of 
the  state  vector  at  t0. 

•  Q  (process  noise  matrix)  -  process  noise  matrix  is 
used  to  define  the  system  model  disturbance  (e.g. 
mean  square  frequency  and  frequency-rate 
fluctuations). 

•  R  (measurement  noise)  -  measurement  noise  matrix 
is  used  to  define  the  measurement  model  based  on 
the  chosen  measurement  errors. 

3.  Initial  conditions  of  the  CSP  are  set:  state  vector  and 
covariance  matrix  which  defines  uncertainty  in  estimate  errors. 

4.  Event  file  (KFTABL.IN)  is  read  containing  one  of  the  four 
events  as  described  above. 

5.  Based  on  the  event  type,  the  next  Kalman  filter  processing  step 
takes  place.  The  normal  procedure  would  be  to  read  a  raw 
measurement,  update  the  covariance  matrix  and  propagate  the 
system  estimates. 

6.  Once  the  measurement  is  read,  an  observation  matrix  H  is 
defined  to  relate  the  measurement  to  the  state  vector. 

7.  The  new  state  vector  elements  are  calculated. 

8.  The  covariance  matrix  diagonal  terms,  phase  and  frequency 
offsets  and  measurement  residuals  are  written  to  an  output  file. 

9.  Steps  5  through  8  are  repeated  until  a  STOP  event  is 
encountered  or  the  final  time  defined  in  the  event  file  is 
reached. 
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BRIEF  DESCRIPTION  OF  INPUT  DATA 


Omega  time  at  each  station  is  derived  from  one  of  three  Hewlett-Packard  CAQI-5061 
Cesium  Beam  Frequency  Standards,  which  are  pan  of  the  AN/FRN-30  Timing-Control  Set.  Each 
frequency  standard  drives  an  Omega  Signal  Format  Generator  (OMSFOG). 

At  present,  all  Omega  transmitting  stations  are  equipped  with  DATUM  model  9390- 
5007  GPS  Time/Frequency  monitors.  USNO  RDAS  is  used  to  gain  access  to  the  input  data  from 
Omega  Stations  Norway  and  Hawaii.  Currently,  at  these  two  stations,  the  RDAS  observes,  records 
and  repons  GPS  measurements  of  UTC  minus  cesium  time  before  it  gets  adjusted  by  OMSFOG, 
as  well  as  the  relative  differences  in  cesium  times.  This  experimental  capability  to  record 
"uncorrected"  data  is  currently  available  only  at  the  Omega  Stations  Norway  and  Hawaii.  A  sample 
of  Precise  Time/Time  Interval  (PTTI)  data  collected  by  USNO  personnel  is  shown  in  Figure  3.2-1, 
an  explanation  of  the  records  is  shown  in  Table  3.2-1.  Each  data  record  field  depicted  in 
Figure  3.2-1  is  described  in  Table  3.2-1  for  Omega  Station  Norway.  Table  3.2-2  contains  record 
field  descriptions  for  PTTI  data  provided  by  Omega  Station  Hawaii.  (Note  that  the  order  of 
measurements  for  primary  and  secondary  cesium  clocks  is  reversed  for  Norway  station:  First, 
secondary  clock  data  is  reported  followed  by  primary  clock  data). 

At  USNO  a  Hewlett-Packard  1000  computer  is  used  to  save  the  PTTI  data  for 
collection  periods  exceeding  one  week.  These  data  were  provided  for  use  in  the  CSP  and  are 
included  on  the  accompanying  diskette  under  filenames  NORWAYl.DAT  and  HAWAII1.DAT. 
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47552  * 

0  * 

.999972641  * 

.000024502  * 

. 015946753 
. 015946998 
.999963657  * 

.000015517  * 


.015937764 
. 015937972 
. 999966471 
. 000018348 
.015940531 
.015940776 
. 399977464 
. 400029324 
.415951522 
.415951756 
99800.3 
-.00000249 


47552.0416667  * 


100  * 

.999972601  * 

.000024463  * 

.  015946627 
.  015946875 
.99996366  * 

.000015521  * 


. 01593768 
. 015937927 
.999966476 
. 00001834 
. 015940503 
. 015940741 
. 399977467 
.40002933 
.415951475 
.415951685 
79700.3 
. 01594707 


Figure  3.2-1  A  Sample  of  PTTI  Data  Collected  by  USNO  RDAS 


*  Direct  input  to  the  CSP 
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Table  3.2-1 

Covariance  Simulation  RDAS  Timing  Measurements 
(Omega  Station  Norway) 


Record  Field 
No. 

Station  Configuration 
Channel  Description 

Description  of  the  Data  Field 

1 

NA 

Modified  Julian  Date  (MJD)  of  the 
data  being  reported  in  this  record* 

2 

NA 

GMT  time  of  the  day* 

3 

Al-Bl 

GPS-CS017* 

4 

A1-B2 

GPS-CS124* 

5 

A1-B3 

GPS-LOR1 

6 

A1-B4 

GPS-LOR2 

7 

A2-B1 

CS486-CS017* 

8 

A2-B2 

CS486-CS124* 

9 

A2-B4 

CS486-LOR1 

10 

A2-B4 

CS486-LOR2 

11 

A3-B1 

OMSFOG-CS017 

12 

A3-B2 

OMSFOG-CS 124 

13 

A3-B3 

OMSFOG-LOR1 

14 

A3-B4 

OMSFOG-LOR2 

15 

NA 

CLOCK-CS047 

16 

NA 

CLOCK-CS124 

17 

NA 

CLOCK-LOR1 

18 

NA 

CLOCK-LOR2 

19 

NA 

CS486-LORM 

20 

! 

NA 

CS486-LORX 

*  Direct  input  to  the  CSP. 
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Table  3.2-2 

Covariance  Simulation  RDAS  Timing  Measurements 
(Omega  Station  Hawaii) 


Record  Field 
No. 

Station  Configuration 
Channel  Description 

Description  of  the  Data  Field 

1 

NA 

Modified  Julian  Date  (MJD)  of  the 
data  being  reported  in  this  record* 

2 

NA 

GMT  time  of  the  day* 

3 

Al-Bl 

GPS-CS554* 

4 

A1-B2 

GPS-CS349* 

5 

A1-B3 

GPS-LOR1 

6 

A1-B4 

GPS-LOR2 

7 

A2-B1 

CS529-CS554* 

8 

A2-B2 

CS529-CS349* 

9 

A2-B4 

CS529-LOR1 

10 

A2-B4 

CS529-LOR2 

11 

A3-B1 

OMSFOG-CS554 

12 

A3-B2 

OMSFOG-CS349 

13 

A3-B3 

OMSFOG-LOR1 

14 

A3-B4 

OMSFOG-LOR2 

15 

NA 

CLOCK-CS554 

16 

NA 

CLOCK-CS349 

17 

NA 

CLOCK- LORI 

18 

NA 

CLOCK-LOR2 

19 

NA 

CS529-LORM 

20 

NA 

CS529-LORX 

*  Direct  input  to  the  CSP. 


Of  the  20  available  data  fields  from  the  USNO  RDAS,  only  the  following  six  data  fields  are 
utilized  by  the  CSP  directly: 

1.  Modified  Julian  Date  (MJD) 

2.  GMT  Time  of  the  Day  (TOD) 

3.  Time  difference  between  GPS  time  and  secondary  cesium  clock 
(GPS-CS017)  for  Norway  Station. 

Time  difference  between  GPS  time  and  primary  cesium  clock  (GPS-CS554)  for  Hawaii 
Station 

4.  Time  difference  between  GPS  time  and  primary  cesium  clock  (GPS-CS124) 
for  Norway  Station. 

Time  difference  between  GPS  time  and  secondary  cesium  clock  (GPS-CS349)  for 
Hawaii  Station 

5.  Time  difference  between  on-line  and  secondary  cesium  clocks  (CS486- 
CS017)  for  Norway  Station. 

Time  difference  between  on-line  and  primary  cesium  clocks  (CS529-CS554)  for  Hawaii 
Station 

6.  Time  difference  between  on-line  and  primary  cesium  clocks  (CS486-CS124) 
for  Norway  Station. 

Time  difference  between  on-line  and  secondary  cesium  clocks  (CS529-CS349)  for 
Hawaii  Station. 

These  are  highlighted  by  an  (*)  in  Table  3.2-1  and  listed  in  the  order  as  they  appear  in 
NORWAY.DAT  file  in  Table  3.2-3.  Table  3.2-4  provides  the  analogous  summary  for 
HAWAII.DAT. 
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Table  3.2-3 

NORWAY.DAT  File  Format 


RECORD 

FIELD 

NO. 

RECORD  CONTENTS 

1 

Modified  Julian  Day 

GMT  Time  of  the  Day 

2 

GPS  -  CS017  (jusec) 

GPS-CS124  (jisec) 

3 

CS486-CS017  (/isec) 

CS486-CS124  (jisec) 

Table  3.2-4 

HAWA13.DAT  File  Format 


RECORD 

FIELD 

NO. 

RECORD  CONTENTS 

1 

Modified  Julian  D 

GMT  Time  of  the  Day 

2 

GPS  -  CS554  (ji sec) 

GPS-CS349  (ftsec) 

3 

CS529-CS554  (/isec) 

CS529-CS349  (/isec) 

3-10 


4.  COVARIANCE  SIMULATION  PROGRAM  OUTPUT 

4.1  SAMPLE  RUNS 

To  demonstrate  the  CSP  performance,  a  discussion  of  three  different  sample  runs  are 
presented  in  this  section.  Each  run  was  defined  by  an  appropriate  update  event  type  that  identifies 
the  measurement  set  and  the  rate  with  which  the  measurement  is  incorporated  into  the  filter.  All 
three  runs  were  executed  on  the  identical  14-day  period  USNO  RDAS  input  data  for  Omega 
Station  Norway  with  a  3-day  scheduled  update  rate.  Table  4.1-1  shows  the  sensors  used  for  the 
sample  runs. 


TABLE  4.1-1 

OMEGA  STATION  NORWAY  SENSOR  UTILIZATION 


SAMPLE  RUN  NO. 

GPS 

TRANSFER 

ON-LINE 

CESIUM 

PRIMARY 

CESIUM 

SECONDARY 

CESIUM 

1 

Off 

On 

On 

On 

2 

On 

On 

Off 

Off 

3 

On 

On 

On 

On 

4.1.1  Sample  Run  1 

Figure  4.1-1  shows  the  plot  of  the  RMS  error  value  of  the  phase  offset  estimate.  This 
plot  was  generated  from  a  data  set  consisting  of: 


•  Time  difference  between  on-line  and  secondary  cesium  measurements 

* 

•  Time  difference  between  on-line  and  primary  cesium  measurements. 


No  GPS  measurement  was  applied  in  this  sample  run.  On  the  third  day  the  same  set  of 
measurements  was  used  to  update  the  uncertainty  of  the  phase  estimate  based  on  the  relative 
measurement  accuracy,  thus  bringing  the  uncertainty  value  down  by  about  0.15  /xsec.  All  three 
cesium  clocks  followed  the  same  "line  of  the  behavior". 


OMEGA  SYNCH  SIMULATION 


T1M90MO 

□  GSAflB  +  0017  o  CS124 

Figure  4.1-1  CS486-CS017,  CS486-CS124  Measurement  Set  (Sample  Run  1) 

4.1.2  Sample  Run  2 

Figure  4.1-2  shows  the  plot  of  the  RMS  error  value  of  the  phase  offset  estimate  for  the 
situation  where  GPS  updates  were  included  for  the  on-line  cesium  standard.  In  this  CSP  run,  it  was 
assumed  that  only  the  on-line  cesium  standard  (CS486)  was  being  compared  to  GPS  time  at 
OMSTA  Norway.  This  sample  run  clearly  demonstrates  the  benefits  of  the  GPS  time  measurement 
incorporation  into  the  system  dynamics.  Even  though  all  three  cesiums  start  out  at  similar  initial 
conditions,  after  the  first  GPS  measurement  update  to  the  on-line  clock,  there  is  a  0.56  jusec 
difference  between  the  uncertainties  of  on-line  and  off-line  cesiums.  As  time  progresses,  this 
uncertainty  appears  to  grow  without  bound  for  the  off-line  cesiums.  It  can  also  be  noted  that 
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incorporating  the  phase  measurements  from  all  three  cesiums  slows  the  rate  of  growth  of  the  phase 
uncertainty.  This  can  be  verified  by  observing  the  value  of  sigmas  in  Figures  4.1*1  and  4.1-2  at  the 
14-day  mark. 

G PS/O MEGA  SYNCH  SIMULATION 

phase  eflVAmucec3-ocr.i  -asm  uk*te) 


CS486 


tiMECDfirs) 

t  C9017 


10  11  12  13  1* 


CS12A 


Figure  4.1-2  GPS-CS486  Measurement  Set  (Sample  Run  2) 


4.1.3  Sample  Run  3 

In  the  final  sample  run,  a  full  station  cesium  setup  was  assumed  to  be  active.  In  other  words, 
all  three  cesium  times  were  used  in  conjunction  with  the  GPS  time  measurement.  The  resulting 
cesium  clock  system  performance  shown  in  Figure  4.1-3  resembles  the  performance  of  a  system 
consisting  of  a  single  on-line  cesium  clock  corrected  by  the  GPS  measurements.  The  combination 
of  the  measurements  used  in  this  run: 


•  Time  difference  between  GPS  and  on-line  cesium 
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•  Time  difference  between  on-line  and  secondary  cesiums 

•  Time  difference  between  on-line  and  primary  cesiums 

As  is  seen  in  Figure  4.1-3,  the  phase  uncertainty  is  bounded  and  well-behaved. 


G PS/OMEGA  SYNCH  SIMULATION 


Figure  4.1-3  GPS-CS486,  CS486-CS017  and  CS486-CS124 
Measurement  Set  (Sample  Run  3) 
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5. 
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APPENDIX  A  FILTER-BAT  FILES 


It  is  imperative  that  the  following  directory  and  their  contents  be  established  in  order  to 
ensure  proper  operation.  The  INSTALL  file  on  the  FILTER  distribution  diskette  will  create  the 
FILTER  directory  and  all  required  files. 


C:\FELTER 


AUTOMENU.COM 

AUTOMENU.MDF 

AUTO.BAT 

FELEPROC.EXE 

GPSSIM.EXE 

PC-WRITE 


FILTER.BAT 

FILTER.MDF 

FDBF.BAT 

FOUT.BAT 

N0RWAY1.DAT 

NORWAY.DAT 

HAWAII1.DAT 

HAWAII.DAT 


C:\LOTUS 

LOTUS.COM 

And  any  application-dependent  files 

I 


FILE  NAMES  AND  FUNCTIONAL  DESCRIPTION 


AUTOMENU.COM 

Executable  menu  application 

AUTOMENU.MDF 

Required  by  AUTOMENU.COM 

AUTO.BAT 

Executes  AUTOMENU.COM 
FILTER.BAT 

Executes  AUTO.BAT  with  FILTER.MDF 
FILTER.MDF 

Contains  definitions  for  FILTER  menu 


A-l 


FDBF.BAT 

Called  by  FILTER.MDF  to  display  a  listing  of  .DAT  and  .BAK 
files  of  the  current  directory 

NORWAY  l.DAT 

USNO  supplied  RDAS  timing  data  for  Omega  Station  Norway 
HAWAII1.DAT 

USNO  supplied  RDAS  timing  data  for  Omega  Station  Hawaii 
NORWAY.DAT 

A  file  containing  converted  Norway  timing  data  (output  of  FELEPROC  program) 
HAWAII.DAT 

A  file  containing  converted  Hawaii  timing  data  (output  of  FILEPROC  program) 
Fileproc  EXE 

Converts  NORWAY  l.DAT  file  into  NORWAY.DAT  file  or  HAWAII1.DAT  file  into 
HAWAII.DAT  file. 

GPSSIM.EXE 

Runs  the  Kalman  filter  simulation  software 
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